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ABSTRACT The action potential of nerve and muscle is produced by voltage-sensitive channels that include a specialized 
device to sense voltage. The voltage sensor depends on the movement of charges in the changing electric field as sug¬ 
gested by Hodgkin and Huxley. Gating currents of the voltage sensor are now known to depend on the movements of posi¬ 
tively charged arginines through the hydrophobic plug of a voltage sensor domain. Transient movements of these 
permanently charged arginines, caused by the change of transmembrane potential V, further drag the S4 segment and 
induce opening/closing of the ion conduction pore by moving the S4-S5 linker. This moving permanent charge induces 
capacitive current flow everywhere. Everything interacts with everything else in the voltage sensor and protein, and so it 
must also happen in its mathematical model. A Poisson-Nernst-Planck (PNP)-steric model of arginines and a mechanical 
model for the S4 segment are combined using energy variational methods in which all densities and movements of charge 
satisfy conservation laws, which are expressed as partial differential equations in space and time. The model computes 
gating current flowing in the baths produced by arginines moving in the voltage sensor. The model also captures the capac¬ 
itive pile up of ions in the vestibules that link the bulk solution to the hydrophobic plug. Our model reproduces the signature 
properties of gating current: 1) equality of ON and OFF charge Q in integrals of gating current, 2) saturating voltage depen¬ 
dence in the Q(charge)-voltage curve, and 3) many (but not all) details of the shape of gating current as a function of 
voltage. Our results agree qualitatively with experiments and can be improved by adding more details of the structure 
and its correlated movements. The proposed continuum model is a promising tool to explore the dynamics and mechanism 
of the voltage sensor. 


INTRODUCTION 

Much of biology depends on the voltage across cell mem¬ 
branes. The voltage across the membrane must be sensed 
before it can be used by proteins. Permanent charges 
move in the strong electric belds within membranes, so car¬ 
riers of sensing charge were proposed as voltage sensors 
even before membrane proteins were known to span lipid 
membranes (1). The movement of permanent charges of 
the voltage sensor is gating current, and the movement is 
the voltage-sensing mechanism. Permanent charge is our 
name for a charge or charge density independent of the 
local electric held (for example, the charge and charge dis¬ 
tribution of Na + but not the charge in a highly polarizable 
anion like Br or the nonuniform charge distribution of 
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H 2 0 in the liquid state with its complex time dependent 
(and perhaps nonlinear) polarization response to the local 
electric held). 

Knowledge of membrane protein structure has allowed 
us to identify and look at the atoms that make up the voltage 
sensor. Protein structures do not include the membrane 
potentials and macroscopic concentrations that power gating 
currents, and therefore, simulations are needed. Atomic-level 
simulations like molecular dynamics (MD) do not provide an 
easy extension from the atomic timescale ~ 10“ 15 s to the bio¬ 
logical timescale of gating currents that starts at ~ 10 6 s 
and reaches ~ 10 2 s. Calculations of gating currents from 
simulations must average the trajectories (lasting ~ 10 ~~ 1 s 
sampled every 10~ 15 s) of ~10 6 atoms, all of which interact 
through the electric held to conserve charge and current 
while conserving mass. It is difficult to enforce continuity 
of current flow in simulations of atomic dynamics because 
simulations compute only local behavior, whereas continuity 
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of current is global, involving current flow far from the atoms 
that control the local behavior. It is impossible to enforce con¬ 
tinuity of current flow in calculations that assume equilib¬ 
rium (zero net flow) under all conditions. 

A hybrid approach is needed, starting with the essential 
knowledge of structure but computing only those parts of 
the structure used by biology to sense voltage. In close- 
packed (“condensed”) systems like the voltage sensor 
or ionic solutions, “everything interacts with everything 
else” because electric fields are long ranged as well as 
exceedingly strong (2). In ionic solutions, ion channels, 
even enzyme active sites, steric interactions that prevent 
the overfilling of space in well-defined protein structures 
are also of great importance because they produce short- 
range correlations (3). 

Closely packed charged systems are well handled math¬ 
ematically by energy variational methods. Energy varia¬ 
tional methods guarantee that all variables satisfy all 
equations (and boundary conditions) at all times and under 
all conditions and are thus always consistent. We use the 
energy variational approach developed in (4) and (5) to 
derive a consistent model of gating charge movement, 
based on the basic features of the structure of crystallized 
voltage-sensitive channels. A schematic of the model is 
shown below. The continuum model we use simulates the 
mechanical dynamics in a single voltage sensor, although 
the experimental data is from many independent voltage 
sensors. Ensemble averages of recordings of individual in¬ 
dependent voltage sensors are equivalent to macroscopic 
continuum modeling in a single voltage sensor if correla¬ 
tions are captured correctly in the model of the single 
voltage sensor. 

MATERIALS AND METHODS 


force (PMF), mainly dominated by the difference of the solvation energy 
in bulk situation and in the hydrophobic plug (7). Note that Na + and Cl - 
(which are the only ions in the bulk solution in this article for simplicity) 
are found only in vestibules and are not allowed into the hydrophobic 
plug in our model. The ends of the two vestibules on each side of the hy¬ 
drophobic plug act as impermeable walls for Na + and Cl - in our model. 
When the voltage is turned on and off, these two walls store/release charge 
(carried by ions) in their electric double layers (EDL) that have many of the 
properties of capacitors. 

In this continuum model, the four arginines (/?,, i = 1, 2, 3, 4) are 
described by their individual density distributions (concentrations) (c,-, 
i — 1, 2, 3, 4), allowing the arginines to interact with Na + and Cl - in 
vestibules. The density (i.e., concentration) distributions represent prob¬ 
ability density functions as shown explicitly in the theory of stochastic 
processes used to derive such equations in (8) using the general methods 
of (9). The important issue here is how well the correlations are captured 
in the continuum model. Some are more likely to be faithfully captured 
in molecular or coarse-grained dynamics simulations (e.g., more or less 
local hard sphere interactions) (10-14) and others in continuum models 
(e.g., correlations induced by far-held boundary conditions like the 
potentials imposed by bath electrodes to maintain a voltage clamp) 
(15-18). 

Here, we treat the S4 itself as a rigid body, so we can capture the basic 
mechanism of a voltage sensor without considering the full details of struc¬ 
ture, which might lead to a three-dimensional model difficult to compute in 
reasonable time. We construct an axisymmetric one-dimensional (ID) 
model with a three-zone geometric configuration illustrated in Fig. 1 b, 
following Fig. 1 a. Zone 1 with z e [0, L R ] is the intracellular vestibule; 
zone 2 with z e [ L R , L R + L\ is the hydrophobic plug; zone 3 with z e 
[L r + L, 2L r + L\ is the extracellular vestibule. Arginines, Na + , and Cl - 
can all reside in zone 1 and 3. Zone 2 only allows the residence of arginines, 
albeit with a severe hydrophobic penalty because of their permanent charge, 
in a region of low dielectric coefficient, hence called hydrophobic. 

Based on Fig. 1 b, the governing ID dimensionless Poisson-Nemst- 
Planck (PNP)-steric equations are expressed below with the detailed nondi- 
mensionalization process shown in Supporting Materials and Methods, 
Section S1. The first one is a Poisson equation that shows how charge cre¬ 
ates potential: 



/ = Na,Cl, 1,2,3,4, (1) 


Theory: Mathematical model 

The reduced mechanical model for a voltage sensor is shown in Fig. 1 a 
with four arginines (R h i— 1, 2, 3, 4), each attached to the S4 helix by iden¬ 
tical springs with the same spring constant K. The electric field will drag 
these four arginines because each arginine carries +1 charge. The charged 
arginines can also move as a group. S4 connects to S3 and S5 at its two ends 
by identical springs with spring constant K S4 /2. 

Once the membrane is depolarized from, for example, —90 mV inside 
negative to +10 mV inside positive, arginines together with S4 will 
be driven toward the extracellular side. A repolarization from +10 to 
—90 mV moves the arginines back to the intracellular side. This movement 
is the basic voltage-sensing mechanism. The movement of S4 triggers the 
opening or closing of the lower gate—consisting mainly of S6 forming 
the ion permeation channel—by a mechanism widely assumed to be me¬ 
chanical, although electrical aspects of the linker motion are likely to be 
involved as well. 

When arginines are driven by an electric field, they are forced to move 
through a hydrophobic plug composed of several nonpolar amino acids 
from SI, S2, to S3 (6). Arginines reside initially in the hydrated lumen of 
the intracellular vestibule. They then move though the hydrophobic plug 
and wind up in the vestibule on the extracellular side. This movement in¬ 
volves dehydration when the arginines move through the hydrophobic 
plug, in which the arginines encounter a barrier in the potential of mean 


where </> is electric potential; is concentration of species i with valence 
qNa ~ 1, gel = -1, qi ~ qarg = 1, i = 1, 2, 3, 4; r = X 2 d /R 2 with 
Xd — >J£ r £§k R TIc^e 1 being the Debye length, and the characteristic length 
(radius of vestibule) R = 1 nm here. A{z) is the channel cross-sectional area 
at position z- For zones 1 and 3, F — 1 by setting NaCl bulk concentration 
Co — 184 mM and £ r — 80 . For zone 2, we assume a hydrophobic environ¬ 
ment with £ r — 8 and therefore F = 0 . 1 . The value of the dielectric constant 
inside the hydrophobic plug (zone 2) is not experimentally available; how¬ 
ever, the computational result is not sensitive to this value based on our 
sensitivity analysis. 

The second equation is the species transport equation based on conserva¬ 
tion laws: 


dcj 

dt 


1 d 
A dz 


(A/,-) = 0, 


i = Na, Cl, 1,2,3,4, 


( 2 ) 


with the content of flux J t expressed below based on the Nemst-Planck 
equation for Na + and Cl - : 


. dCi <9</A 

Ji = —D, ——|- Ciqr— I, / = Na, Cl, z in zone 1 and 3, 
dz dz / 


(3) 
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FIGURE 1 (a) Geometric configuration of 

gating pore in this model, including the attach¬ 
ments of arginines to the S4 segment. ( b ) Following 
(a), an axisymmetric three-zone domain shape is 
designated in r-z coordinate for the current ID 
model. Here, the diameter of the hydrophobic 
plug is 0.3 nm (arginine’s diameter); L — 0.7 nm; 
L r = 1.5 nm; and the radius of the vestibule is 
R — 1 nm. BC means boundary condition. To see 
this figure in color, go online. 


and for four arginines c„ i = 1, 2, 3 and 4 based on the Nemst-Planck equa¬ 
tion with steric effect and some imposed potentials: 



where Z), is the diffusion coefficient for species i. 

The first and second terms in Eqs. 3 and 4 describe diffusion and electro¬ 
migration, respectively. The third terms in Eq. 4 are external potential terms 
with Vi, i— 1,2, 3, and 4 being the constraint potential for the four arginines 
Ci to S4, represented here by a spring connecting each arginine c, to S4, as 
shown in Fig. 1 a. Governing equations Eqs. 1, 2, 3, and 4 were derived by 
energy variational methods, which is further shown in Supporting Materials 
and Methods, Section S3. 

The elastic system is described by 

V/M = K(z — (z,- + Zs 4 (f))) 2 , (5) 

where K is the spring constant, n is the fixed anchoring position of the 
spring for each arginine c z on S4, and Z s4 (t) is the center-of-mass z position 


of S4 by treating S4 as a rigid body. Here, we set zi — 0.6, Z 2 — 0.2, 
Z 3 = —0.2, and z 4 = —0.6 using structural information that gives the argi¬ 
nine anchoring interval on S4 as 0.4 nm. Z S4 (t) follows the motion of equa¬ 
tion based on the spring-mass system: 


d 2 Z<:A dZ?A 

m S4 ~E3-1" ^54 , + Ks4 (Zs4 — Zs4, o) 


dt 2 


dt 


4 

— K( z i,CM — (Zi + Z 54 )), 
i — 1 


( 6 ) 


where m S4 , b S4 , and K S4 are the mass, damping coefficient, and restraining 
spring constants for S4. Zs4,o is the resting position of Zs 4 (t). Here, Zi,cM is 
the center of mass for the set of arginines c h which can be calculated by 


r 

Jo 


•L+2LR 


z i,CM r L+2LR 


A(z)zc;dz 

W ,i = 1, 2, 3, 4. (7) 


foA(z)c,dz 


We assume that the spring-mass system for S4 is overdamped, which 
means the inertia term in Eq. 6 can be neglected. 

The energy barrier V b in Eq. 4 is nonzero only in zone 2, which mainly 
represents the difference in solvation energy, chiefly characterized by the 
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difference of dielectric constants, in the hydrophobic plug and bulk solu¬ 
tion. The structure of the energy barrier is actually very complicated. 
Here, we simply assume a hump shape for PMF (see more in Supporting 
Materials and Methods, Section S2), although we will seek greater realism 
in later work. 

The last term in Eq. 4 is the steric term that accounts for steric interaction 
among arginines (5,19). Here, we set g — 0.5, a reasonable value. Though 
there is actually no experimental measurement available for g, the compu¬ 
tation results have been verified to be insensitive to its value. 

Here, we assume quasisteady state for Na + and Cl - , which means 
dci/dt = 0, i = Na, Cl, in Eq. 2, and the reasons are elaborated in Support¬ 
ing Materials and Methods, Section S4. The formulation of boundary and 
interface conditions is also shown in Supporting Materials and Methods, 
Section S5. 

Besides the main input parameter V, which is the applied voltage bias 
(corresponding to the command potential in voltage-clamp experiments), 
other parameters like D, (i — 1, 2, 3, 4), K, K S4 , and b S4 are also required. 
Results are especially sensitive to the values of K, K$ 4 , and bs 4 . We 
have tried and found D, = 50; i — 1, 2, 3, and 4; K — 3; K S4 = 3; 
and bs 4 =1.5 provide the best fit to the experimental Q(charge)-voltage 
(QV) curve reported in (20). Some additional explanation on fitting these 
parameter values is described in Supporting Materials and Methods, 
Section S6. 

Usually, the electric current in the ion channel is treated simply as the 
flux of charge and is uniform in the z direction when steady in time. This 
is not so in this nonsteady dynamic situation because the storing and 
releasing of charge in vestibules is involved. Here, the flux of charge 
at the middle of hydrophobic plug, z = L R + LI 2, was computed to es¬ 
timate the experimentally observed gating current. However, it is actu¬ 
ally impossible (so far) to experimentally measure the current at the 
middle of the hydrophobic plug. In experiments, the voltage-clamp tech¬ 
nique is used, and on/off gating current through the membrane is 
measured, which should be equal to the flux of charge at z — 0 in this 
framework, as shown in Fig. 1 b. The flux of charges at any z position 
7(z, t ) can be related to the flux of charges at z — 0, 1 ( 0 , t ), simply by 
charge conservation: 

J<2net(z,t) = 1(0, t) ~I(z,t), (8) 

where 

z 

Q ne ,(z,t) = f A(q)'y^q i c I dE, ( 9 ) 

o a " 1 

and flux of charges at any z position I(z, t) is defined by 

/(V) = AWJXv). (10) 

all i 

We identify d/dtQ ne[ (z,t) as the displacement current and denote it as 
Idisp ( z , t) because Eq. 8 is equivalent to Ampere’s law in Maxwell’s equa¬ 
tions, and d / dt(Q nel (z, t)) is exactly the displacement current in Ampere’s 
law. The proof is elaborated on in Supporting Materials and Methods, Sec¬ 
tion S7. A general discussion about displacement current can be found in 
(21-23), which does not involve assumptions concerning the dielectric co¬ 
efficient e r or polarization properties of matter at all. Hence, Eq. 8 can be 
simply rewritten as 

I,o,{z,t) = I(z,t) +I disp (z,t) = 1(0 ,t), ( 11 ) 

where we define the sum of displacement current and flux of charges as 
the total current I tot (z, t). The z distribution of the total current should be 


uniform by Kirchhoff’s law, and we verify this by computations shown 
in the section under heading Flux of Charges at Different Locations. 
Note the ionic current I(z, t) changes a great deal with location. The 
displacement current Idi sp (z, t) varies a great deal with location. The total 
current, the sum I tot (z, t), does not vary at all with location, although of 
course it varies a great deal with time. For example, calculations of cur¬ 
rent in the baths (which are not reported here) would show only ionic 
current in the time range considered here, but it would equal the total 
current that flows anywhere in our ID model of the voltage sensor 
domain. 

We are also interested in observing the net charge at vestibules. Consider, 
for example, the net charge at the intracellular vestibule, Q ne t(^R , t). The net 
charge consists of arginine charges and their countercharges formed by the 
EDL of ionic solution in that location. Electroneutrality is approximate but 
will not be exact there. Flux of charge, displacement current, and net charge 
at vestibules will be discussed further in the section under heading Flux of 
Charges at Different Locations. 

To evaluate the current theoretical model, it is important to compare 
our computational results with experimental measurements (20) in 
the curves of gating current and amount of gating charge moved 
versus applied voltage (I(current)-voltage [IV] and QV curves). 
To construct the QV curve, we calculate Q\ — f Q R A(z)J2f=iC,dz, Q 2 — 
f£* +L A(z)Y^ i= \Cidz, Q 3 — A(z)Y^ = i c idz, which are the amounts 

of arginine found in zone 1, 2, and 3, respectively. Usually Q 2 ~ 0 is due 
to the energy barrier Vj, in zone 2. Arginines tend to jump across zone 2 
when driven from zone 1 to zone 3 as the voltage V is turned on. The 
number of arginines that move and settle at zone 3 depends on the 
magnitude of V. Besides IV and QV curves, the time course of the move¬ 
ment of arginines and S4, Zi,cM(t) and Zs 4 (t), is important to report here 
because recording these movements in experiments is becoming feasible 
nowadays by optical methods. Many qualitative models accounting for 
the movement of S4 and conformation change of the voltage sensor 
have been proposed. Readers are referred to review articles (24,25) for 
more details. 


Numerical method 

Equations 1, 2, 3, and 4 are first discretized in space by high-order multi¬ 
block Chebyshev pseudospectral methods and then integrated in time under 
the framework of method of lines. The details of the numerical method are 
referred to Supporting Materials and Methods, Section S8. 


RESULTS AND DISCUSSION 

Here, numerical results based on the mathematical model 
described above were calculated and compared with exper¬ 
imental measurements (20). Our ID continuum model has 
advantages and disadvantages. The lack of three-dimen¬ 
sional structural detail means that some details of the gating 
current and charge cannot be reproduced. It should be noted, 
however, that to reproduce those, one needs more than just 
static structural detail. One must also know how the struc¬ 
tures (particularly their permanent and polarization charge) 
move and change after a command potential is applied in the 
experimental ionic conditions. The ID model has advan¬ 
tages because it computes the actual experimental results 
on the actual experimental timescale in realistic ionic solu¬ 
tions and with far-held boundary conditions actually used in 
voltage-clamp experiments. It also conserves total current, 
as we will demonstrate later. Conservation of current needs 
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to be there and verified in theories and simulations because 
it is a universal property of the Maxwell equations (21-23). 

QV curve 

When the membrane and voltage sensor is held at a large in¬ 
side negative potential (e.g., hyperpolarized to —90 mV), S4 
is in a resting potential position, and all arginines stay in the 
intracellular vestibule. When the potential is made more 
positive (e.g., depolarized to +10 mV), S4 is in the active 
potential position, and all arginines are at the extracellular 
vestibule. 

The voltage dependence of the charge (arginines) trans¬ 
ferred from intracellular vestibule to extracellular vestibule 
is characterized as a QV curve in experimental papers, and 
it is sigmoidal in shape (20). Fig. 2 a shows that our 
computed QV curve—the dependence of Qi on V —is in 
very good agreement with the experiment (20). This good 
agreement comes from the fact that our resultant QV curve 
is also a sigmoidal curve, and, most important of all, the 
slope of QV curve can be tuned, mainly by the adjustment 
of K , K s4 , and b s4 , to agree with experiment. Not many theo¬ 
retical models can achieve this agreement, especially for the 
slope. Models in (15,16) show good agreement with exper¬ 
iments, whereas a mismatch of slope was reported in 
(17,18). The voltage dependence of activation has been 
considered a crucial property of the sodium conductance 
since it was defined (1). Fig. 2 b shows the steady-state 
distributions of Na + , Cl - , and arginines in the inside nega¬ 
tive, hyperpolarized situation (V = —90 mV). As we can 
see, all the arginines stay in the intracellular vestibule, 
and none of the arginines move to the extracellular vestibule 
(23 « 0). 

Fig. 2 c shows the situation at V = —48 mV, which is the 
midpoint of the QV curve. As we can see, each vestibule has 
distributions of c t (i = 1, 2, 3, and 4), resulting in half of the 
arginines staying in it (23 = 2). The center-of-mass position 
for each arginine, presented later in Fig. 6, shows that R1 
and R2 are in the extracellular vestibule, and R3 and R4 


are in the intracellular vestibule. There are almost no argi¬ 
nines in zone 2 (hydrophobic plug) because of the energy 
barrier in it. Note that this represents an average because 
in a single molecule interpretation, half of the sensors will 
be with all R’s inside and the other half with all R’s outside. 
The midpoint of —48 mV from (20) requires the resting po¬ 
sition of S4, ^54,0- to be biased from L R + 0.5 L to Zs4,0 — 
L r + 0.5 L + 1.591 nm; otherwise, the midpoint would be 
0 mV. Fig. 2 cl shows the situation at full depolarization 
(V = —8 mV), at which time all arginines move to the extra¬ 
cellular vestibule (23 ~ 4) in the fully depolarized, acti¬ 
vated state. 

Gating current 

Fig. 3 shows the time course of gating currents, observed 
as flux of charge at the middle of hydrophobic plug 
I(L r + LI 2, t) because of the movement of arginines when 
the membrane depolarization is large and when the depolar¬ 
ization is small. In the case of large depolarization, V rises 
from —90 mV at t = 10 to —8 mV and drops back to 
—90 mV at t = 150 (Fig. 3 a). The time course of gating cur¬ 
rent and contributions of individual arginines are shown in 
Fig. 3 b. As expected, the rising order of each current 
component follows the moving order of Rl, R2, R3, and 
R4 when depolarized and that order is reversed when repo¬ 
larized. The area under the gating current is the amount of 
charge moved. Because arginines move forward and back¬ 
ward in this depolarization/repolarization scenario, the areas 
under the ON current and the OFF current are same. The 
areas are equal for each component of current as well. 
The equality of area is an important signature of gating 
current that contrasts markedly with the properties of ionic 
current (26,27). In the case of small depolarization (V rises 
from —90 to —50 mV at t— 10 to and drops back to —90 mV 
at t = 150, Fig. 3 c), the time course of gating current and its 
four components contributed by each arginine for this situ¬ 
ation is shown in Fig. 3 cl. Under this small depolarization, 
not all arginines move past the middle of the hydrophobic 



0 12 3 

z (nm) 


0 12 3 

z (nm) 


FIGURE 2 (a) QV curve and comparison 

with (20). Steady-state distributions for Na + , Cl - , 
and arginines are shown at ( b ) V = —90 mV, 
(c) V — —48 mV, and ( d) V — —8 mV. Note that 
the experimental data in (20) were scaled to 4e. To 
see this figure in color, go online. 
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FIGURE 3 (a) The time course of V rising from —90 to —8 mV at t = 10 holds on until t = 150 and then drops back to —90 mV. (b) The time course of 

gating current, I(L R + LI 2, /), and its components corresponding to (a) are shown, (c) The time course of V rising from —90 to —50 mV at t = 10 holds on until 
t = 150 and then drops back to —90 mV. (d) The time course of gating current. I(L R + LI 2, /). and its components corresponding to (c) are shown. To see this 
figure in color, go online. 


plug because of the weaker driving force in the small depo¬ 
larization compared with the large depolarization case. This 
can be inferred because the areas under each component 
current are different (Fig. 3 d). 

The gating currents can be better understood by looking 
at a sequence of snapshots showing the spatial distribution 
of electric potential, species concentration, and electric 
current. The distributions at several times are shown in 
Fig. 4 a for the case of sudden change in command voltage 
to a more positive value and a large depolarization, and the 
distributions are shown in Fig. 4 b for the case of a small 
depolarization. The electric potential profiles at t = 13 
and t = 148 show that the profile of electric potential 
changes as arginines move from left to right even though 
the voltage is maintained constant across the sensor. Slight 
bulges in electric potential profile exist wherever arginines 
are dense. This can be easily explained by understanding 
the effect of Eq. 1 on a concave spatial distribution of elec¬ 
tric potential. 

In Fig. 4, the total current defined in Eq. 11, though 
changing with time, is always constant in z at all times, satis¬ 
fying Kirchhoff’s law (i.e., conservation of current). At 
t = 13, when gating current is substantial, as seen from 
t = 13 in Fig. 3, b and <7, we can visualize the z distributions 
of flux of charges I(z, t). displacement of current l,jj sp (z, t), 
and total current I tot {z, t) individually in Fig. 4. 

Flux of charges at different locations 

Flux of charges I(z, t), together with displacement current 
Idispiz , t) and total current t), depicted in Fig. 4, deserve 

more discussion here. Though I(z, t), Idispiz , f), and I tot (z, t) are 


well defined inEqs. 8,9, 10, and 11, the actual computation of 
them takes an indirect path because of the assumption of qua¬ 
sisteady state for Na + and Cl - in Eq. 2. The details are pre¬ 
sented in Supporting Materials and Methods, Section S9. 
The computed total current I tot (z, t) does indeed satisfy 
Kirchhoff’s law by its uniformity in z. This verification is 
shown in Fig. 4 at several times, and we have checked that 
this is in fact true at any time. 

In the bottom rows of Fig. 4 at t = 13, we observe that 
I(z , t) is generally nonuniform in z and is accompanied by 
congestion/decongestion of arginines in between. However, 
I(z, t) is almost uniform in zone 2 (hydrophobic plug), which 
means almost no congestion/decongestion of arginines oc¬ 
curs there, and therefore, there is no contribution to the 
displacement current d/dtQ nel {z,t) from zone 2. This is 
because arginines can hardly reside in zone 2 because of 
the energy barrier in it. 

Several things are worth noting in the time courses of 
I{Lr + L/2, r) and 7(0, t) (equal to uniformly distributed 
I tot as depicted by Eq. 11) illustrated in Fig. 5 a under the 
case of large depolarization. First, I(L R + L/2, t) is notice¬ 
ably larger than 7(0, t) in the ON period. This is because their 
difference, exactly the displacement current I disp , is always 
negative at zone 2 when depolarized because arginines 
are leaving zone 1 and make d/dtQ net < 0 for zone 2. 
It is expected that the area under the time course of 
I(L r + L/2, r) would be very close to 4e, as verified 
by the time courses of Q 3 in Fig. 5 b. We use 7(0, t ) to esti¬ 
mate the experimentally measured voltage-clamp current, 
whereas the counterpart area of experimentally measurable 
7(0, t ) would be less than 4e because of its smaller magni¬ 
tude compared with 7(L S + L/2,f). This may partly explain 
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FIGURE 4 (a) The top row shows dimensionless 
species concentration distributions at t = 0, 13 
(right after depolarization), and 148 (right before 
repolarization) for the case of large depolarization 
with V from —90 mV at t — 10 to —8 mV and drop¬ 
ping back to —90 mV at t — 150. The middle 
row shows concurrent electric potential profiles. 
The bottom row shows concurrent electric current 
profiles with the components of flux of charge, 
displacement current, and total current. ( b ) The 
same as (a) is shown except with V depolarized 
from —90 to —50 mV. To see this figure in color, 
go online. 


the experimental observations that at most 13e (25,28,29), 
instead of 16e, are moved during full depolarization in 
four voltage sensors (for a single ion channel) based on 
computing the area under voltage-clamp gating current. 
Therefore, flux of charge at any location of zone 2, though 
impossible to measure in experiments so far, will give us 
the amount of arginines moved during depolarization 
more reliably than the measurable 7(0, t). 

Second, we see in Fig. 5 a with magnification in its inset 
plot that, as in experiments, 7(0, t), but not7(L s + L/2,t), has 
contaminating leading spikes in ON and OFF parts of the 
current. These spikes are capacitive currents from solution 
EDL of vestibules caused by the sudden rising and dropping 
of command potential. These spikes need to be removed in 
voltage-clamp experiments to get rid of the contribution 
from vestibule solution EDL (and membrane) to the trans¬ 
port of gating charges (arginines) when computing the 
area under gating current. The technical details of removing 
these spikes are shown in Supporting Materials and 


Methods, Section S10, and more details about spikes can 
be found in Supporting Materials and Methods, Section SI 1. 

Third, in Fig. 5 b, as arginines move from one vestibule 
to another, the concentrations of Na + and Cl - also corre¬ 
spondingly change with time at the vestibules. They form 
countercharges through EDL and balance arginine charges 
at vestibules. However, these EDL changes only maintain 
an approximate, not exact, charge balance, as shown in 
Fig. 5 b. The violation of electroneutrality causes the 
displacement current, which is not negligible. This further 
causes the underestimate of arginines that move when the 
voltage sensor is depolarized if the estimate is made by 
measuring the area under 7(0, t). 

As in the previous section, we used flux of charges at the 
middle of the hydrophobic plug, I(L R + 7/2, t), instead of 
experimentally measurable 7(0, t) to represent the gating 
current in discussions. We may as well name I(L R + 
772, t) as the arginine current to avoid the confusion with 
the actual gating current 7(0, t ) here. This arginine current 
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FIGURE 5 (a) The time courses of I(L R + L/2,t), 

1 (0, t), and despiked 1 (0, t) for the case of large depo¬ 
larization with V rising from — 90 to — 8 mVat f = 10 , 
holding on till t = 150 , and then dropping back 
to —90 mV. The inset plot is a magnification 
of the ON current to visualize the difference of 
1(0, t) and despiked 1(0, t ) more clearly, (h) The 
time courses of Q\ . Q\. J f} 11 (c Ra — Cci)dz, and 
fu+L L ( CNa ~ c ci)dz are under the same depolariza¬ 
tion scenario as (a). To see this figure in color, 
go online. 


leaves out its associated displacement current Idi sp (L R + 
LI 2, 1) and serves to represent gating current better for two 
reasons: 

1) The area under the time course of I(L R + L/2, t) gives us 
the amount of arginines moved during depolarization 
more faithfully than 1(0, t). The fluxes of charge for 
each arginine shown in Fig. 3, b and d carry important 
information about how each arginine is moved by the 
electric field that will be further illustrated in Fig. 6. 
All these will not be easy to display and comprehend if 
we use 7(0, t) instead. 

2) Using 1(0, t) as a definition of gating current would 
require a decontamination by removing the leading 
spikes, which is computationally costly. Removing 
spikes would especially pose a heavy numerical burden 
when doing parameter fitting in which numerous 
repeated computations are done. 

Time course of arginine and S4 translocation 

Fig. 6 shows the time course of Q (amount of arginines 
moved to extracellular vestibule, equal to <2 3 here) and cen- 


ter-of-mass trajectories of individual arginines (Zi,cM, 7=1, 
2, 3, and 4) and S4 segment (Z s4 ). Fig. 6, a and b show the 
case of large depolarization, and Fig. 6, c and d show the 
case of small depolarization. 

In the case of large depolarization (Fig. 6 b), the argi¬ 
nines and S4 z positions quickly reach individual steady 
states, with almost all arginines transferred to the extracel¬ 
lular vestibule as previously shown in Fig. 4 a. Therefore, 
Q is close to its saturated value 4 as shown in Fig. 6 a. Ar¬ 
ginines and S4 move back to the intracellular vestibule once 
the voltage drops back to —90 mV. From Fig. 6 b, the for- 
ward-moving order of arginines is Rl, R2, R3, and R4, and 
the backward-moving order is the opposite R4, R3, R2, 
and Rl with agreement with the structure. This agreement 
might look trivial in molecular dynamics simulations but 
is not a trivial check here because this model describes 
arginines not by particles, as in molecular dynamics, but 
by concentrations. Note that an incorrect order and pace 
of the movement of arginines would cause disagreement 
with experiments in the shape of IV curve as well. S4 is 
initially farthest to the right but lags behind Rl and R2 dur¬ 
ing movement in depolarization, as shown in Fig. 6 b. This 
is certainly because S4 is finally relaxed to an almost 
unforced situation close to its resting position Z.S4,0 during 



FIGURE 6 ( a ) and (c) are the time courses of the 
amount of arginines moved to the extracellular 
vestibule. ( b ) and ( d) are center-of-mass trajec¬ 
tories of individual arginines and S4. ( a ) and ( b ) 
are the case of large depolarization with V rising 
from —90 to —8 mV at t = 10, holding on till 
t — 150, and then dropping back to —90 mV. 
(c) and ( d) are the case of small depolarization 
with V rising from —90 to —50 mV at t — 10, hold¬ 
ing on till t — 150, and then dropping back to 
—90 mV. To see this figure in color, go online. 
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this large depolarization. We can further calculate the 
displacements of each arginine and S4 during this full- 
saturating depolarization and find Az\,cm ~ Az 2 ,cm ~ 
Az%,cm ~ 1-93 nm, Az 4 ,cm = 1-76 nm, and AZ S4 = 
1.51 nm. Besides almost the same displacements for Rl, 
R2, and R3, their average moving velocities are also very 
close to each other. This seems to suggest a synchronized 
movement among Rl, R2, and R3 that we have not imposed 
on the arginines in our model. Also, we can see the move¬ 
ments of arginines contribute significantly to the movement 
of the S4 segment. This can be seen from the steady-state z 
position of S4 derived from Eq. 6, 


_ K -A K S4 

Zsi — — -—rr; } y (Zj,CM - Z/J + -p -TTiT^ 4 . 0 


Km + 4K 

1 


i=l 

4 


K S4 + 4 K 


Zs4,0 + Zi ,i 


i,CM 


(ID 


Experimental estimates of S4 displacement during full de¬ 
polarization range from 2 to 20 A (24,30), depending on the 
model of the voltage sensor and its motion, including the 
transporter model, the helical screw, and the paddle model 
(24). Our AZ S4 =1.51 nm here is large and seems to agree 
better with experimental estimates requiring large displace¬ 
ments, such as the paddle model. In contrast, the helical 
screw model, which is supported by most of the recent 
data, is known to have shorter displacements. A plausible 
explanation for our overestimate of AZ S4 is that our ID model 
uses a straight line perpendicular to the hydrophobic-plug 
path for the movement of the arginines. In reality, the S4 
segment is significantly tilted with respect to the membrane, 
and the arginines follow a spiral along the helix. Therefore, if 
the S4 segment rotates and changes its tilt during activation, 
the total vertical translation needed to cross the hydrophobic 
plug is significantly reduced, as was shown by Vargas et al. 


(31). The value obtained in (31) was between 0.7 and 1 nm 
when comparing the displacement perpendicular to the mem¬ 
brane of the open-relaxed state crystal structure of Kv 1.2 (32) 
and the closed structure that has been derived by consensus 
from experimental measurements (31). 

In the case of small depolarization, the driving force is 
weaker than in a large-saturating depolarization, so their 
Z positions do not have a chance to reach steady states as 
they do during a full-saturating depolarization. Rather, 
in a small depolarization, the motion of the arginines and 
S4 are aborted. They return to the intracellular vestibule 
because the depolarization drops (i.e., decreases in magni¬ 
tude, and the membrane potential becomes more negative) 
before arginines and S4 have a chance to reach their 
steady-state positions. This detailed atomic interpretation 
likely overreaches the resolution of our model. At the sin¬ 
gle-sensor level, we do not expect partial movements; 
instead, some sensors will have moved all the way and 
others not at all, but the distribution of sensors in the two 
extreme positions should follow what we predict with this 
model, which is an ensemble average. We look forward to 
measurements of movements of probes that mimic arginine 
in its environment that require improvements in the resolu¬ 
tion and structural realism of our model. 

Fig. 6 c illustrates these aborted motions. Q reaches 1.57 
at most, which should be 2 instead if steady-state was 
reached as it is if time is long enough. See the steady-state 
behavior shown in the QV curve of Fig. 2 a. Fig. 6 cl shows 
that the S4 segment is initially farthest to the right, lags 
behind Rl during movement, and is almost caught up by 
R2. The maximal displacements of arginines and S4 calcu¬ 
lated from Fig. 6 d are Az\,cm — 1-36 nm, Az 2 ,cm = 
0.966 nm, Azi,cm = 0.459 nm, Az 4 ,cm = 0.316 nm, and 
AZ 4cm = 0.616 nm. The significant difference between 
Azi,cm, Az 2 .cm, Az X cm , and Az 4 , C m may imply that Rl 
and R2 have jumped across the hydrophobic plug and 
entered the extracellular vestibule, whereas R3 and R4 


a v (mV) b 



FIGURE 7 (a) The time courses of subtracted 
gating current, despiked 1 (0, t), with the voltage 
rising from —90 to V mV at t = 10, holds on till 
t = 150, and then drops back to —90 mV, where 
V = —62, —50, ..., —8 mV (b) r 2 versus V 
compared with experiment (20) is shown. To see 
this figure in color, go online. 
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still stay at the intracellular vestibule during this small 
depolarization. This is consistent with the observation 
from individual gating-current components of arginines in 

Fig. 3 d. 

Family of gating currents for a range of voltages 

Though we prefer 1(Lr + L/2, t) to 7(0, t) for representing 
gating current as explained in the section under heading 
Flux of Charges at Different Locations, we here use the 
actual gating current, despiked 7(0, t), to compare with 
experiment (20). Fig. 7 a shows the time courses of a sub¬ 
tracted gating current (despiked 7(0, t)) for a range of volt¬ 
ages V ranging from —62 to —8 mV. The area under 
gating current, for both ON and OFF parts, increases 
with V because more arginines are transferred to the extra¬ 
cellular vestibule as V increases. The shapes of this family 
of gating currents agree well with experiment (20) in both 
magnitude and time course. 

We can characterize the time course by fitting the decay 
part of a subtracted gating current by ae~'' T ' + be~'' T2 , 
t i < t 2 as generally done in experiments (20) in which t, 
is the fast time constant and t 2 is the slow time constant. 
Usually, the movement of arginines is dominated by t 2 . 
Flere, t 2 was calculated from simulation and compared 
with experiment (20) as shown in Fig. 7 b. Because in our 
computation the time is in arbitrary units, we have scaled 
the time to have the maximal t 2 to fit with its counterpart 
in experiment (20). Overall, the trend of t 2 versus V in our 
result, though not the whole curve, agrees well with exper¬ 
iment (20). To the left of the maximal point in Fig. 7 b, 
simulation results fit rather well with the experiment 
compared with the values to the right of the maximal point, 
at which it overestimates t 2 compared with the experiment. 
This overestimate is consistent with the observation that the 
amount of transferred charges Q saturates slightly faster in 
experimental data than in this simulation as V increases 
(see QV curve of Fig. 2 a). This phenomenon is related 


to the cooperativity of movement among arginines, which 
will be further discussed below. 

Effect of voltage pulse duration 

Fig. 8 shows the effect of voltage pulse duration with Fig. 8 a 
for the case of small depolarization and Fig. 8 b for the case of 
large depolarization. The magnitude and time span of sub¬ 
tracted gating current (despiked 7(0, t )) are changed by pulse 
duration in both cases, but the shape will asymptotically 
approach the same curve as pulse duration increases, no 
matter the size of the depolarization. This behavior occurs 
because it takes time for the command pulse to drive the ar¬ 
ginines toward the extracellular vestibule. If the pulse dura¬ 
tion is long enough, the time course of Q will approach its 
steady state for large depolarization as in Fig. 6 a. Small de¬ 
polarization takes a longer time to reach its steady state, as 
demonstrated in Fig. 6 c. The shapes of gating currents in 
Fig. 8 compare favorably with experiment (20) in which 
the OFF subtracted gating currents for short pulses have 
very fast decays, whereas for long pulses, the OFF subtracted 
gating currents have larger rising amplitude and slower decay 
because of a larger amount of arginines moved. 

CONCLUSIONS 

Previous work with molecular and coarse-grained simula¬ 
tions have captured some interactions, but they have not 
yet reproduced the time course and voltage dependence of 
macroscopic gating currents (10-14), and previous contin¬ 
uum models have captured only the steady-state properties 
of charge movement (15-18). 

This ID continuum mechanical model of the voltage 
sensor tries to capture the essential structural details of the 
movement of mass and charge that are necessary to repro¬ 
duce the basic features of experimentally recorded gating 
currents. After finding appropriate parameters, we find 
that the general kinetic and steady-state properties are 




FIGURE 8 Subtracted gating currents, despiked 
7(0, t ), showing the effect of voltage pulse duration. 
(a) V increases from —90 to —35 mV at t — 10 and 
drops back to —90 mV at various times. ( b ) V in¬ 
creases from —90 to 0 mV at t — 10 and drops 
back to —90 mV at various times. To see this figure 
in color, go online. 
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well represented by the simulations. The good agreement of 
our numerical results with salient features of gating current 
measured experimentally would be impossible by simply 
tuning of parameters if our model had not captured the 
essence of physics for the voltage sensor. The continuum 
approach seems to be a good model of voltage sensors, pro¬ 
vided that it 1) takes into account all interactions crucial to 
the movement of gating charges and S4; 2) computes their 
correlations consistently, so all variables satisfy all equa¬ 
tions under all conditions with one set of parameters; 
and 3) satisfies conservation of current. This last point 
gave us a new insight: what is measured experimentally 
does not correspond to the transfer of the arginines because 
the total current, containing a displacement current, is 
smaller than the arginine current. It should be noted, how¬ 
ever, that the total energy provided by the voltage clamp 
is qV, where q is the time integral of the measured gating 
current and Vis the applied voltage. This is the total energy 
that explains the correspondence of charge per channel with 
the charge estimated by the limiting slope method ( 33 - 35 ). 

We have simplified the profile of the energy barrier in the 
hydrophobic plug because the PMF in that region, and its 
variation with potential and conditions, is unknown. There 
is plenty of detailed information on the amino acid side 
chains in the plug and how each one of them changes the ki¬ 
netics and steady-state properties of gating charge move¬ 
ment (6). Therefore, the next step is to model the details 
of interactions of the moving arginines with the wall of 
the hydrophobic plug and the contributions from other sur¬ 
rounding charged protein components. Some of the effects 
to be included are as follows: 

1) Steric and dielectric interactions of the arginines that this 
model does not include. These include the interaction of 
arginines with negative charges of the S2 and S3 seg¬ 
ments and the negative phospholipids as well as the hy¬ 
drophobic residues in the plug. These interactions may 
be responsible for the simultaneous movement of two 
to three arginines across the plug, which is an experi¬ 
mental result that this model does not reproduce ( 36 , 37 ). 

2) Time dependence of the plug energy barrier V*. Once the 
first arginine enters the hydrophobic plug by carrying 
some water with it, this partial wetting of the hydropho¬ 
bic plug will lower Vt, chiefly consisting of solvation en¬ 
ergy, and enable the next arginine to enter the plug with 
less difficulty. This might explain the cooperativity of 
movement among arginines when they jump through 
the plug. The addition of details in the plug may also pro¬ 
duce intermediate states that have been measured exper¬ 
imentally. In this situation, arginines may transiently 
dwell within the plug. 

3) A very strong electric field might affect the hydration 
equilibrium of the hydrophobic plug and would lower 
its hydration energy barrier as well ( 38 ). This cooper¬ 
ativity of movement may help explain the quick satura¬ 


tion in the upper right branch of the QV curve (and 
smaller r 2 ). It may also explain the experimentally 
observed translocation of two to three arginines simulta¬ 
neously ( 36 , 37 ). 

The power of this mathematical modeling is precisely the 
implementation of interactions and the various effects in a 
consistent manner. Implementing the various effects listed 
above is likely to lead to a better prediction of the currents 
and to the design of experiments to further test and extend 
the model. 

Further work must address the mechanism of coupling 
between the voltage sensor movements and the conduction 
pore. For example, the spring constant of the two sides of 
S4 have been made equal, which does not take into account 
the structural reality that one side has a linker to S3, 
whereas the other links to the pore opening. It seems likely 
that the classical mechanical models of coupling will need 
to be extended to include coupling through the electrical 
field. The charges involved are large. The distances are 
small, so the changes in electric forces that accompany 
movements of charged mass (and flows of displacement 
current) are likely to be large and important. It is possible 
that the voltage sensor modifies the stability of a funda¬ 
mentally stochastically unstable, nearly bistable, conduc¬ 
tion current (of single channels) by triggering sudden 
transitions from closed to open state in a controlled 
process reminiscent of Coulomb blockade in a noisy 
environment ( 39 ). 
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1. Non-dimensionalization 


We non-dimensionalize all physical quantities as follows, 


c~ =f,4> = 
c 0 


JJ = — 

k B T/e’ k B T ’ 


S = t = 


R 2 /d x 


Di = 


Bi g _ 9ij I _ Jj 

D x ’ k B T/c 0 ’ Jl c 0 D x /R’ 


I = 


I 


ec 0 D x R 


where q is concentration of species /, with /=Na + , Cl”, 1, 2, 3, and 4. Each is 
scaled by c 0 which is the bulk concentration of NaCl in the 
intracellular/extracellular domains. Here c 0 is set to be 184 mM, equal on both 


sides, so that the Debye length A D = 


e r E 0 k B T 


is lnm when the relative 


»/ Coe 2 

permittivity e r = 80. <fi is the electric potential scaled by k B T /e with k B 
being the Boltzmann constant; T the temperature; e the elementary charge. All 
relevant external potentials U are scaled by k B T. All sizes s are scaled by R, which 
is the radius of vestibule as shown in Fig. 1(b). f?=lnm here. The time t is scaled 
by R 2 /D x , with D x being a diffusion coefficient that can be adjusted later to be 
consistent with the time spans of on/off currents measured in experiments 
(caused by the movement of arginines). The diffusion coefficient of species / is 
scaled by D x . The coupling constant of PNP-steric model based on 
combining rules of Lennard Jones, representing the strength of steric interaction 
between species / and j, is scaled by k B T/c 0 [1,2]. For simplicity, we assume 

= [ 0 ’ ffor^ll/i' 5 —y ' = Note that here we only consider steric 

interaction among arginines. We think they are a crucial source of correlated 
structural change and motion (of mass and charge). The consideration of steric 
effect among arginines is justified by the fact that arginines are generally 
crowded in hydrophobic plug and vestibules. The flux density of species /, Ji, is 
scaled by c 0 D x /R, and therefore the electric current / is scaled by ec 0 D x R. For 
simplicity of notation, we will drop ~ for all dimensionless quantities shown in 
all equations. 


2. Shape of potential of mean force (PMF) in the hydrophobic plug 

Here, we simply assume a hump shape for PMF in the hydrophobic plug as, 


Vb = V b max (tanh(5(z — L R )) — tanh(5(z — L — L R )) — l), when z is in zone 2, 
Vb = 0, when z is in zone 1 and 3, 


(SI) 


with V b max set to be 5 for a good agreement with experimental measurements. 


l 



Theoretically, if we set V b max too large, the gating current would be slow and 
perhaps small because it would be very difficult for arginines to move across this 
barrier. The double tanh functions are designed to smooth the otherwise 
top-hat-shape barrier profile, which is not good for numerical differentiation 
because of its awkward infinite slopes. This smoothing is simply based on the 
belief that the energy barrier in a protein structure does not have a jump. In 
future work, it would be wise to compute the PMF from a specific model of charge 
distribution (both permanent and polarization) constructed from a combination 
of structural data and molecular dynamics simulations, if feasible. 


3. Governing equations derivation from energy variation methods 

Governing equations Eqs. (1-4) were derived by energy variational methods 
based on the following energy (in dimensional form): 


E = L k B T Yjall i c ilog c i - I V01 2 + Yuall i Ri e c i0 + £ 


arginines 


(v i + v b )c i + 


j arginines 


9ij 

■ ■ —- c r ■ 
IJ 2 1 J. 


dV, 


(S.2) 


where the first term is entropy; second and third terms are electrostatic energy; 
the fourth term is the constraint and barrier potential for arginines; the last term 
is the steric energy term, based on Lennard-Jones potential [1,3]. The Poisson 
equation Eq. (1) is derived from the variation of energy with respect to electric 
potential 


SE 



and species flux densities in Eqs. (3,4) are derived by 


SE 

~sF t ‘ ,l ~ 


o, 


kuT 




where /ij is the chemical potential of species /. 


4. Quasi-steadiness assumption for Na + and Cl 

9 C ■ 

Here we assume quasi-steady state for Na + and CP, which means = 0, i = 

Na, Cl. The steady state assumption here is justified by the fact that the diffusion 
coefficients of Na + and Cl in vestibules are much larger than the diffusion 
coefficient of arginine based on the very narrow time span of the leading spike of 
gating current measured in experiments. The spike comes from the linear 
capacitive current of vestibule when the command potential suddenly rises or 
drops. This quasi-steady state assumption is essential for the success of our 


2 



calculations. Otherwise using realistic diffusion coefficients for Na + and Cl-would 
render Eqs. (1-4) too stiff to integrate in time. The spike contaminating the gating 
current is removed in experiments by a simple technique called P/n leak 
subtraction (see Section 11; n typically is 4). P/n leak subtraction is also used to 
subtract the linear capacity current of all the membranes in the real system that 
are not included in our model. How to do leak subtraction computationally will 
be discussed in Section 10. 


5. Formulation of boundary conditions 

Types of boundary conditions are illustrated in Fig. 1(b). Note the no-flux 
boundary conditions specified in Fig. 1(b). One prevents Na + and Cl from 
entering the hydrophobic plug (zone 2) with low dielectric coefficient. The other 
boundary condition constrains S4 motion and so prevents the arginines from 
leaving the vestibules into intracellular/extracellular domains. 

Boundary and interface conditions for electric potential <fi are 

dd) dd> 

0 ( 0 ) = v, (P(Lr) = <pa + R ), r(LfiM(Lfi)-^(LK) = r(L+M(L+)-^(L+), 

4kl r + n = 0(L fl + l + ), t(l* + l-)a(l r + + n = + 


L + )A(L r + L + ) g (Lr + L + ), 0(2L R + L) = 0. (S.3) 

These are Dirichlet boundary conditions at both ends and continuity of electric 
potential and displacement at the interfaces between zones. Boundary and 
interface conditions for arginine are 

Ji(. 0> t) = Ji(/L R + L,t) = 0, Ci(L R ,t ) = Q (L r , t), A(L R )J i (L R , t) = 

A(L R )Ji(L R ,t ), Ci(L r + L~,t) = Ci(L r + L + ,t), A(L R + L~)J/L R + L~,t ) = 

A(L r +l + )Ji{l R +L + ,t), i = 1,2,3,4, (S.4) 

where no-flux boundary conditions are placed at both ends of the gating pore, 
consisting of vestibules and hydrophobic plug, to prevent arginines and S4 from 
entering intracellular/extracellular domains. The others are continuity of 
concentration and flux at interfaces between zones. Boundary conditions for Na + 
and CF are 


c W a(0,0 = c a (0,t) = c Na (2L R + L,t) = c cl (2L R + L,t) = 1, 

= Jci(L R ,t) = Jnci(Lr + L,t ) = Jci(L r + L, t) = 0, (S.5) 


where Dirichlet boundary conditions are placed at both ends of the gating pore to 
describe the concentrations for Na + and Cl as the bulk concentration. No-flux 
boundary conditions at both ends of hydrophobic plug describe the 


3 



impermeability of Na + and Cl into hydrophobic plug. 


6. Parameters fitting 

We have tried and found D,=50, /= 1,2,3,4, K= 3, K S4 = 3, b S 4 = 1.5 provide the 
best fit to the important experiments reported in [4]. Several things are to be 
noted about the parameter values specified above: (1) there is no experimental 
measurement of diffusion coefficient of arginine inside vestibule and plug 
available that we can use for simulation. Imprecise setting of the values of these 
diffusion coefficients only affects the scale of time in I-V curve, but not its shape. 
That is why we set time coordinate to be in an arbitrary unit later in results, and 
here we only focus on comparing the shape of IV curves with experiments in [4]. 
(2) K, Ks 4 , and bs 4 were particularly determined by fitting with QV curve in 
experiment [4]. The QV curve is very sensitive to K and Ks 4 , and many efforts have 
been taken to achieve proper values for them. The method of fitting is done by 
trial and error. Choosing incorrect K and Ks 4 would end up serious mismatch of 
QV curve with experiment [4] as demonstrated by the case of K= 3 and Ks 4 = 12 in 
Fig. 1 here. The choice of K= 3 and Ks 4 = 3 fits experiment [4] best and is adopted 
for the rest of simulations. 

4 

3 

O 2 

1 

0 

Figure 1. Simulated QV curves under different K and Ks 4 compared with 
experimental counterpart from [4]. Note that the experimental data in [4] was 
scaled to 4e. 

7. Derivation of Ampere's law in Maxwell's equations by Poisson equation 
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and species transport equation 

Eq. (8) is consistent with Ampere's law in Maxwell's equations: 

Vx © = £ » £ -37 + '- ( S - 6 ) 

or equivalently, 

V'(^ r f+/) = 0, (S.7) 

where E is the electric field and / is flux density of charge (current density). Eq. 
(S.7) tells us that the total current is conserved everywhere and it consists of flux 

of charges / and displacement current e 0 £ r —. Eq. (S.7) can be derived from the 

Poisson equation and species transport equation like Eq. (1) and Eq. (2). Starting 
from Poisson equation in dimensional form: 

-V ■ (s 0 £ r V(p) = p + It qieci, (S.8) 
or equivalently 

V ■ ( £ 0 £ r E ) = p + Ej qteci. (S.9) 

Taking time derivative of Eq. (S.9), 
v-( Vr f)=Z i <? i ^, CS.io) 

and using species transport equation based on mass conservation, 

$ + v-/, = 0, (S.11) 

then 

V • ( £ o £ r ft. ) = & 1 ‘ e S = - V ■ Si Wh = -V ■ /. (S.12) 

which becomes exactly Eq. (S.7) by defining 

/ = Si <7t e /r (S-13) 

A more general treatment that does not involve assumptions about £ r can be 
found in [5-7]. 

Casting Eq. (S.7) into the present ID framework by integrating it in space and 
applying the divergence theorem, we have 

£ 0 £ r A(z ) d -^ + I(z, t) = £ 0 £ r A(0) d -^ + 1(0, t). (S.14) 

Comparing with Eq. (11), 

s o s rA(z) ^ - £ 0 £ r A( 0) ^ = I dlsp (z, t), (S.15) 

which justifies the naming of displacement current in Eq. (11). 
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8. Numerical method 

High-order multi block Chebyshev pseudospectral methods are used here 
to discretize Eqs. (1-4) in space [8]. The resultant semi discrete system is then a 
set of coupled ordinary differential equations in time and algebraic equations (an 
ODAE system) [9]. The ordinary differential equations are chiefly from Eq. (2), 
and algebraic equations are chiefly from Eq. (1) and boundary/interface 
conditions Eqs. (S.3-S.5). This system is further integrated in time by an ODAE 
solver (ODE15S in MATLAB (The Math Works, Natick, MA) [10,11]) together with 
appropriate initial condition. ODE15S is a variable order variable step (VSVO) 
solver, which is highly efficient in time integration because it adjusts the time step 
and order of integration. High order pseudospectral methods generally provide 
excellent spatial accuracy with economically practicable resolutions. A 
combination of these two techniques makes the whole computation very efficient. 
This is particularly important here, since numerous computations have to be 
tried during the tuning of parameters. Efficiency will be vital in future 
calculations comparing theory and experiment in a wide variety of mutants and 
experimental conditions. 

9. Computation of flux of charge, displacement current and total current 

According to definition in Eq. (10), flux of charges at the middle of gating 
pore, I(L r + L/2, t), and both ends of gating pore, 7(0, t) and I(2L R + L, t), 
should be computed by 

1 (I'R T ~, t) A ^ L r + Z larginines RiJi (j^R T (S.16) 

7(0, t) = A{ 0) Y.i=Na,ci qdii 0,0, (S.17) 

I(2L r +L,t) = A(2L r + L) Z i=Na , cl qJi(2L R + L, t). (S.18) 

Except I[L R +^,t^, 7(0, t) and I(2L R + L,t) are trivially zero due to the 

d c ■ 

implement of quasi-steadiness = 0, i = Na, Cl, in vestibules, which causes 

J Na and J C i to be uniform in vestibules by Eq. (2), and further become zero by 
the no-flux boundary conditions for Na + and Cl - at the bottom of vestibules as 
described in Eq. (S.5). We have to alternatively reconstruct 7(0, t) and 
I(2L R + L,t) by charge conservation of Na + and Cl - , 

K°,£) = A(z)Y.,i„.aq i c l dz l (S.19) 

1(2L r + L, t) = -^S L l ^ a (z) Zm,a Widz. (S.20) 
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After obtaining 7(0, t) and I(2L R + L, t), we can further reconstruct the flux of 
charges 7(z, t) at zone 1 and zone 3 by (8) and (9), 

I(z,t) = l(0,t) ~Y t fo M^ZaiuChCidz, z e [0,L R ], (S.21) 

7(z, t) = I(2L r + L, t) + ^ J z 2ii?+i A(z) S aa t q lCl dz, z E [L R + L, 2 L R + L\. (S.22) 
Flux of charge at zone 2 is simply 

I(z, t) t1(z) ^arginines RiJi^> 0» ^ ^ [Lp, L R + L], (S.23) 

since Na + and Cl" are not allowed to enter zone 2, the hydrophobic plug. 

10. Removing spike in total current 

In voltage-clamp experiments, subtracting this linear capacitive component 
and removing the spike from gating current is done by 'leak subtraction', in 
various forms, e.g., P/4 (see details in Section 11) In reality, this linear capacitive 
current that is subtracted in this procedure comes from both the lipid bilayer 
membrane in parallel with the gating pore. Here, we only considered the 
capacitive current from solution EDL of vestibule inside the gating pore and 
ignored the membrane capacitive current because we simply use Dirichlet 
boundary conditions for <fi at both ends of the gating pore in Eq. (S.3). Actually, 
capacitive current of the membrane in parallel with the gating pore would be 
much larger than vestibule capacitive current. Following the idea of the 
experiment [4], we calculated 7(0, t) with V rising from -150 mV to -140 mV at 
t=10, and dropping back to -150 mV at t=150. We chose from -150 mV to -140 
mV because essentially none of the arginines move across the hydrophobic plug 
in this hyperpolarized region. The voltage step quickly charges and discharges 
solution EDL in vestibules, and the computed time course of 7(0, t) is just two 
spikes at on and off of the command potential. Subtracting this hyperpolarized 
7(0, t), multiplied by a proportion factor (due to the linearity of capacitive 
current), from its original counterpart will then remove the spikes, and the 
unspiked 7(0, t) is shown in Fig. 5(a). In preliminary calculations with the model, 
when the command voltage pulse rises faster, the early spike becomes larger and 
is still visible even after subtraction, suggesting that is the origin of the early 
transient gating current in experiments [12-14]. 

11. Removing linear capacitive current to obtain gating current in 
experiments 

Our computations have limited fidelity at short times because of time step 
limitations in integrating stiff systems. The spike artifacts are one example, 
described previously. Experimental measurements [12,15] of the fast transient 


7 



gating current are fascinating and our calculations will be extended to explore 
more of them in future study by using greater resolution in time. 

A more general consideration is the subtraction procedure used in 
experiments to isolate gating current from currents arising from other sources. 
Channels and their voltage sensors are embedded in lipid membranes, therefore 
they are 'in parallel' with large capacitive currents of the lipid bilayer. The lipid 
membrane has a large capacitance ( C lipid = 8x 10 -7 farads/cm 2 ) that has 
nothing to do with the current produced by charge movement in the voltage 
sensor. Fortunately, the capacitance C Upid is a nearly ideal circuit element and the 
current to charge it is entirely a displacement current accurately described by 
i cap = Qtptd dV/dt with a single constant C lipid . V is the voltage across the lipid 
capacitor. Note that i cap does not include any current or flux of charge carried 
across the lipid. 

In experimental measurements, i cap is always present. Experimental 
measurements always mix the displacement currents of lipid membrane and 
voltage sensor. Lipid membrane current usually dominates the measurement of 
gating currents in native preparations and remains large in systems mutated to 
have unnaturally large numbers of voltage sensors. 

A procedure to remove the lipid membrane current is needed if the gating 
current of the voltage sensor is to be measured. The procedure introduced by [16] 
has been used ever since in the improved P/4 version developed by [17] 
reviewed and discussed in [18]. Also, see another approach in [19] and [20]. 
Schneider and Chandler's procedure [19] estimates the so-called linear current 
i x = C x dV/dt in conditions in which the voltage sensor and C x behave as ideal 
circuit elements. The voltage sensor might then have a component linear in 
potential. An ideal capacitor has a capacitance C x independent of voltage, time, 
current, or ionic composition. The Schneider procedure then subtracts that linear 
current i x —plus any linear component of voltage sensor current—from the total 
current measured in conditions in which the voltage sensor does not behave as 
an ideal capacitor. The leftover estimates the nonlinear properties of the charge 
movement in the voltage sensor. That is to say, the leftover estimates the charge 
movement of the voltage sensor that is not proportional to the size of the voltage 
step used in the measurement. The leftover is called gating current here and in 
experimental papers. 

The gating current reported in experiments [16] can miss a component of the 
displacement current of the voltage sensor, if it uses the linear subtraction to 
estimate i x . These procedures can remove more than the current through the 
lipid membrane capacitor i cap . 
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Clearly, some of the current produced by movements of the arginines in the 
voltage sensor will be a linear displacement current, a linear component of gating 
current and it would not be present in the reported gating current determined by 
some linear subtraction procedures. In particular, if the arginine system is 
present at the 'control' potential contributing a current linear in potential, this 
problem would occur. Of course, if the arginine system is immobilized and 
inactivated at the control potential and so contributes no current flow under that 
condition, this problem would not occur. 

Other systems may contribute to the linear displacement current as well, for 
example, i) all sorts of experimental and instrumentation artifacts and ii) 
displacement current in the conduction channel itself. The conduction channel of 
field effect transistors produces a large displacement current often characterized 
as a capacitance that involves diffusion and is described by drift diffusion 
equations quite similar to the PNP equations of the open conduction channel. 

Most systems have substantial motions that are linear in voltage (even if the 
system is labeled 'nonlinear'). The linear term is present in most systems, just as 
it is present in most Taylor expansions of nonlinear functions. 

The linear component that can be missed in experiments, and removed in 
these calculations, may have functional and structural significance. The voltage 
sensor works by sensing voltage, for example, by producing a motion of arginines. 
That motion—the response of the voltage sensor in this model—includes a linear 
component. The signal passed to the conduction channel, to control gating, is 
likely to include or depend on the linear component of sensor function. Confusion 
will result if a significant linear component exists and is ignored when a model is 
created that links the voltage sensor to the gating process of the conduction 
channel. Direct measurements of the movement of arginines (e.g., with optical 
methods) are likely to include the linear component and so should not agree with 
experimental measurements of gating current or with the currents reported here 
if the linear component exists and is significant in size. 

If the P/4 procedure subtracts a charge movement in a control system in 
which the arginines do not move at all (because they are immobilized and 
inactivated, in that sense), then the resulting estimate of gating current will 
contain a component linear in voltage. Thus, the interpretation of the corrected 
record depends on the details of immobilization and inactivation, topics that are 
beyond the scope of this paper and our present work. 
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c~ =f,4> = 
c 0 


JJ = — 

k B T/e’ k B T ’ 


S = t = 


R 2 /d x 


Di = 


Bi g _ 9ij I _ Jj 

D x ’ k B T/c 0 ’ Jl c 0 D x /R’ 


I = 


I 


ec 0 D x R 


where q is concentration of species /, with /=Na + , Cl”, 1, 2, 3, and 4. Each is 
scaled by c 0 which is the bulk concentration of NaCl in the 
intracellular/extracellular domains. Here c 0 is set to be 184 mM, equal on both 


sides, so that the Debye length A D = 


e r E 0 k B T 


is lnm when the relative 


»/ Coe 2 

permittivity e r = 80. <fi is the electric potential scaled by k B T /e with k B 
being the Boltzmann constant; T the temperature; e the elementary charge. All 
relevant external potentials U are scaled by k B T. All sizes s are scaled by R, which 
is the radius of vestibule as shown in Fig. 1(b). f?=lnm here. The time t is scaled 
by R 2 /D x , with D x being a diffusion coefficient that can be adjusted later to be 
consistent with the time spans of on/off currents measured in experiments 
(caused by the movement of arginines). The diffusion coefficient of species / is 
scaled by D x . The coupling constant of PNP-steric model based on 
combining rules of Lennard Jones, representing the strength of steric interaction 
between species / and j, is scaled by k B T/c 0 [1,2]. For simplicity, we assume 

= [ 0 ’ ffor^ll/i' 5 —y ' = Note that here we only consider steric 

interaction among arginines. We think they are a crucial source of correlated 
structural change and motion (of mass and charge). The consideration of steric 
effect among arginines is justified by the fact that arginines are generally 
crowded in hydrophobic plug and vestibules. The flux density of species /, Ji, is 
scaled by c 0 D x /R, and therefore the electric current / is scaled by ec 0 D x R. For 
simplicity of notation, we will drop ~ for all dimensionless quantities shown in 
all equations. 


2. Shape of potential of mean force (PMF) in the hydrophobic plug 

Here, we simply assume a hump shape for PMF in the hydrophobic plug as, 


Vb = V b max (tanh(5(z — L R )) — tanh(5(z — L — L R )) — l), when z is in zone 2, 
Vb = 0, when z is in zone 1 and 3, 


(SI) 


with V b max set to be 5 for a good agreement with experimental measurements. 
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Theoretically, if we set V b max too large, the gating current would be slow and 
perhaps small because it would be very difficult for arginines to move across this 
barrier. The double tanh functions are designed to smooth the otherwise 
top-hat-shape barrier profile, which is not good for numerical differentiation 
because of its awkward infinite slopes. This smoothing is simply based on the 
belief that the energy barrier in a protein structure does not have a jump. In 
future work, it would be wise to compute the PMF from a specific model of charge 
distribution (both permanent and polarization) constructed from a combination 
of structural data and molecular dynamics simulations, if feasible. 


3. Governing equations derivation from energy variation methods 

Governing equations Eqs. (1-4) were derived by energy variational methods 
based on the following energy (in dimensional form): 


E = L k B T Yjall i c ilog c i - I V01 2 + Yuall i Ri e c i0 + £ 


arginines 


(v i + v b )c i + 


j arginines 


9ij 

■ ■ —- c r ■ 
IJ 2 1 J. 


dV, 


(S.2) 


where the first term is entropy; second and third terms are electrostatic energy; 
the fourth term is the constraint and barrier potential for arginines; the last term 
is the steric energy term, based on Lennard-Jones potential [1,3]. The Poisson 
equation Eq. (1) is derived from the variation of energy with respect to electric 
potential 


SE 



and species flux densities in Eqs. (3,4) are derived by 


SE 

~sF t ‘ ,l ~ 


o, 


kuT 




where /ij is the chemical potential of species /. 


4. Quasi-steadiness assumption for Na + and Cl 

9 C ■ 

Here we assume quasi-steady state for Na + and CP, which means = 0, i = 

Na, Cl. The steady state assumption here is justified by the fact that the diffusion 
coefficients of Na + and Cl in vestibules are much larger than the diffusion 
coefficient of arginine based on the very narrow time span of the leading spike of 
gating current measured in experiments. The spike comes from the linear 
capacitive current of vestibule when the command potential suddenly rises or 
drops. This quasi-steady state assumption is essential for the success of our 


2 



calculations. Otherwise using realistic diffusion coefficients for Na + and Cl-would 
render Eqs. (1-4) too stiff to integrate in time. The spike contaminating the gating 
current is removed in experiments by a simple technique called P/n leak 
subtraction (see Section 11; n typically is 4). P/n leak subtraction is also used to 
subtract the linear capacity current of all the membranes in the real system that 
are not included in our model. How to do leak subtraction computationally will 
be discussed in Section 10. 


5. Formulation of boundary conditions 

Types of boundary conditions are illustrated in Fig. 1(b). Note the no-flux 
boundary conditions specified in Fig. 1(b). One prevents Na + and Cl from 
entering the hydrophobic plug (zone 2) with low dielectric coefficient. The other 
boundary condition constrains S4 motion and so prevents the arginines from 
leaving the vestibules into intracellular/extracellular domains. 

Boundary and interface conditions for electric potential <fi are 

dd) dd> 

0 ( 0 ) = v, (P(Lr) = <pa + R ), r(LfiM(Lfi)-^(LK) = r(L+M(L+)-^(L+), 

4kl r + n = 0(L fl + l + ), t(l* + l-)a(l r + + n = + 


L + )A(L r + L + ) g (Lr + L + ), 0(2L R + L) = 0. (S.3) 

These are Dirichlet boundary conditions at both ends and continuity of electric 
potential and displacement at the interfaces between zones. Boundary and 
interface conditions for arginine are 

Ji(. 0> t) = Ji(/L R + L,t) = 0, Ci(L R ,t ) = Q (L r , t), A(L R )J i (L R , t) = 

A(L R )Ji(L R ,t ), Ci(L r + L~,t) = Ci(L r + L + ,t), A(L R + L~)J/L R + L~,t ) = 

A(L r +l + )Ji{l R +L + ,t), i = 1,2,3,4, (S.4) 

where no-flux boundary conditions are placed at both ends of the gating pore, 
consisting of vestibules and hydrophobic plug, to prevent arginines and S4 from 
entering intracellular/extracellular domains. The others are continuity of 
concentration and flux at interfaces between zones. Boundary conditions for Na + 
and CF are 


c W a(0,0 = c a (0,t) = c Na (2L R + L,t) = c cl (2L R + L,t) = 1, 

= Jci(L R ,t) = Jnci(Lr + L,t ) = Jci(L r + L, t) = 0, (S.5) 


where Dirichlet boundary conditions are placed at both ends of the gating pore to 
describe the concentrations for Na + and Cl as the bulk concentration. No-flux 
boundary conditions at both ends of hydrophobic plug describe the 
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impermeability of Na + and Cl into hydrophobic plug. 


6. Parameters fitting 

We have tried and found D,=50, /= 1,2,3,4, K= 3, K S4 = 3, b S 4 = 1.5 provide the 
best fit to the important experiments reported in [4]. Several things are to be 
noted about the parameter values specified above: (1) there is no experimental 
measurement of diffusion coefficient of arginine inside vestibule and plug 
available that we can use for simulation. Imprecise setting of the values of these 
diffusion coefficients only affects the scale of time in I-V curve, but not its shape. 
That is why we set time coordinate to be in an arbitrary unit later in results, and 
here we only focus on comparing the shape of IV curves with experiments in [4]. 
(2) K, Ks 4 , and bs 4 were particularly determined by fitting with QV curve in 
experiment [4]. The QV curve is very sensitive to K and Ks 4 , and many efforts have 
been taken to achieve proper values for them. The method of fitting is done by 
trial and error. Choosing incorrect K and Ks 4 would end up serious mismatch of 
QV curve with experiment [4] as demonstrated by the case of K= 3 and Ks 4 = 12 in 
Fig. 1 here. The choice of K= 3 and Ks 4 = 3 fits experiment [4] best and is adopted 
for the rest of simulations. 

4 

3 

O 2 

1 

0 

Figure 1. Simulated QV curves under different K and Ks 4 compared with 
experimental counterpart from [4]. Note that the experimental data in [4] was 
scaled to 4e. 

7. Derivation of Ampere's law in Maxwell's equations by Poisson equation 


experiment 

simulation, K=3, K S4 =3 
simulation, K=3, K S4 =12 
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and species transport equation 

Eq. (8) is consistent with Ampere's law in Maxwell's equations: 

Vx © = £ » £ -37 + '- ( S - 6 ) 

or equivalently, 

V'(^ r f+/) = 0, (S.7) 

where E is the electric field and / is flux density of charge (current density). Eq. 
(S.7) tells us that the total current is conserved everywhere and it consists of flux 

of charges / and displacement current e 0 £ r —. Eq. (S.7) can be derived from the 

Poisson equation and species transport equation like Eq. (1) and Eq. (2). Starting 
from Poisson equation in dimensional form: 

-V ■ (s 0 £ r V(p) = p + It qieci, (S.8) 
or equivalently 

V ■ ( £ 0 £ r E ) = p + Ej qteci. (S.9) 

Taking time derivative of Eq. (S.9), 
v-( Vr f)=Z i <? i ^, CS.io) 

and using species transport equation based on mass conservation, 

$ + v-/, = 0, (S.11) 

then 

V • ( £ o £ r ft. ) = & 1 ‘ e S = - V ■ Si Wh = -V ■ /. (S.12) 

which becomes exactly Eq. (S.7) by defining 

/ = Si <7t e /r (S-13) 

A more general treatment that does not involve assumptions about £ r can be 
found in [5-7]. 

Casting Eq. (S.7) into the present ID framework by integrating it in space and 
applying the divergence theorem, we have 

£ 0 £ r A(z ) d -^ + I(z, t) = £ 0 £ r A(0) d -^ + 1(0, t). (S.14) 

Comparing with Eq. (11), 

s o s rA(z) ^ - £ 0 £ r A( 0) ^ = I dlsp (z, t), (S.15) 

which justifies the naming of displacement current in Eq. (11). 
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8. Numerical method 

High-order multi block Chebyshev pseudospectral methods are used here 
to discretize Eqs. (1-4) in space [8]. The resultant semi discrete system is then a 
set of coupled ordinary differential equations in time and algebraic equations (an 
ODAE system) [9]. The ordinary differential equations are chiefly from Eq. (2), 
and algebraic equations are chiefly from Eq. (1) and boundary/interface 
conditions Eqs. (S.3-S.5). This system is further integrated in time by an ODAE 
solver (ODE15S in MATLAB (The Math Works, Natick, MA) [10,11]) together with 
appropriate initial condition. ODE15S is a variable order variable step (VSVO) 
solver, which is highly efficient in time integration because it adjusts the time step 
and order of integration. High order pseudospectral methods generally provide 
excellent spatial accuracy with economically practicable resolutions. A 
combination of these two techniques makes the whole computation very efficient. 
This is particularly important here, since numerous computations have to be 
tried during the tuning of parameters. Efficiency will be vital in future 
calculations comparing theory and experiment in a wide variety of mutants and 
experimental conditions. 

9. Computation of flux of charge, displacement current and total current 

According to definition in Eq. (10), flux of charges at the middle of gating 
pore, I(L r + L/2, t), and both ends of gating pore, 7(0, t) and I(2L R + L, t), 
should be computed by 

1 (I'R T ~, t) A ^ L r + Z larginines RiJi (j^R T (S.16) 

7(0, t) = A{ 0) Y.i=Na,ci qdii 0,0, (S.17) 

I(2L r +L,t) = A(2L r + L) Z i=Na , cl qJi(2L R + L, t). (S.18) 

Except I[L R +^,t^, 7(0, t) and I(2L R + L,t) are trivially zero due to the 

d c ■ 

implement of quasi-steadiness = 0, i = Na, Cl, in vestibules, which causes 

J Na and J C i to be uniform in vestibules by Eq. (2), and further become zero by 
the no-flux boundary conditions for Na + and Cl - at the bottom of vestibules as 
described in Eq. (S.5). We have to alternatively reconstruct 7(0, t) and 
I(2L R + L,t) by charge conservation of Na + and Cl - , 

K°,£) = A(z)Y.,i„.aq i c l dz l (S.19) 

1(2L r + L, t) = -^S L l ^ a (z) Zm,a Widz. (S.20) 
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After obtaining 7(0, t) and I(2L R + L, t), we can further reconstruct the flux of 
charges 7(z, t) at zone 1 and zone 3 by (8) and (9), 

I(z,t) = l(0,t) ~Y t fo M^ZaiuChCidz, z e [0,L R ], (S.21) 

7(z, t) = I(2L r + L, t) + ^ J z 2ii?+i A(z) S aa t q lCl dz, z E [L R + L, 2 L R + L\. (S.22) 
Flux of charge at zone 2 is simply 

I(z, t) t1(z) ^arginines RiJi^> 0» ^ ^ [Lp, L R + L], (S.23) 

since Na + and Cl" are not allowed to enter zone 2, the hydrophobic plug. 

10. Removing spike in total current 

In voltage-clamp experiments, subtracting this linear capacitive component 
and removing the spike from gating current is done by 'leak subtraction', in 
various forms, e.g., P/4 (see details in Section 11) In reality, this linear capacitive 
current that is subtracted in this procedure comes from both the lipid bilayer 
membrane in parallel with the gating pore. Here, we only considered the 
capacitive current from solution EDL of vestibule inside the gating pore and 
ignored the membrane capacitive current because we simply use Dirichlet 
boundary conditions for <fi at both ends of the gating pore in Eq. (S.3). Actually, 
capacitive current of the membrane in parallel with the gating pore would be 
much larger than vestibule capacitive current. Following the idea of the 
experiment [4], we calculated 7(0, t) with V rising from -150 mV to -140 mV at 
t=10, and dropping back to -150 mV at t=150. We chose from -150 mV to -140 
mV because essentially none of the arginines move across the hydrophobic plug 
in this hyperpolarized region. The voltage step quickly charges and discharges 
solution EDL in vestibules, and the computed time course of 7(0, t) is just two 
spikes at on and off of the command potential. Subtracting this hyperpolarized 
7(0, t), multiplied by a proportion factor (due to the linearity of capacitive 
current), from its original counterpart will then remove the spikes, and the 
unspiked 7(0, t) is shown in Fig. 5(a). In preliminary calculations with the model, 
when the command voltage pulse rises faster, the early spike becomes larger and 
is still visible even after subtraction, suggesting that is the origin of the early 
transient gating current in experiments [12-14]. 

11. Removing linear capacitive current to obtain gating current in 
experiments 

Our computations have limited fidelity at short times because of time step 
limitations in integrating stiff systems. The spike artifacts are one example, 
described previously. Experimental measurements [12,15] of the fast transient 
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gating current are fascinating and our calculations will be extended to explore 
more of them in future study by using greater resolution in time. 

A more general consideration is the subtraction procedure used in 
experiments to isolate gating current from currents arising from other sources. 
Channels and their voltage sensors are embedded in lipid membranes, therefore 
they are 'in parallel' with large capacitive currents of the lipid bilayer. The lipid 
membrane has a large capacitance ( C lipid = 8x 10 -7 farads/cm 2 ) that has 
nothing to do with the current produced by charge movement in the voltage 
sensor. Fortunately, the capacitance C Upid is a nearly ideal circuit element and the 
current to charge it is entirely a displacement current accurately described by 
i cap = Qtptd dV/dt with a single constant C lipid . V is the voltage across the lipid 
capacitor. Note that i cap does not include any current or flux of charge carried 
across the lipid. 

In experimental measurements, i cap is always present. Experimental 
measurements always mix the displacement currents of lipid membrane and 
voltage sensor. Lipid membrane current usually dominates the measurement of 
gating currents in native preparations and remains large in systems mutated to 
have unnaturally large numbers of voltage sensors. 

A procedure to remove the lipid membrane current is needed if the gating 
current of the voltage sensor is to be measured. The procedure introduced by [16] 
has been used ever since in the improved P/4 version developed by [17] 
reviewed and discussed in [18]. Also, see another approach in [19] and [20]. 
Schneider and Chandler's procedure [19] estimates the so-called linear current 
i x = C x dV/dt in conditions in which the voltage sensor and C x behave as ideal 
circuit elements. The voltage sensor might then have a component linear in 
potential. An ideal capacitor has a capacitance C x independent of voltage, time, 
current, or ionic composition. The Schneider procedure then subtracts that linear 
current i x —plus any linear component of voltage sensor current—from the total 
current measured in conditions in which the voltage sensor does not behave as 
an ideal capacitor. The leftover estimates the nonlinear properties of the charge 
movement in the voltage sensor. That is to say, the leftover estimates the charge 
movement of the voltage sensor that is not proportional to the size of the voltage 
step used in the measurement. The leftover is called gating current here and in 
experimental papers. 

The gating current reported in experiments [16] can miss a component of the 
displacement current of the voltage sensor, if it uses the linear subtraction to 
estimate i x . These procedures can remove more than the current through the 
lipid membrane capacitor i cap . 
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Clearly, some of the current produced by movements of the arginines in the 
voltage sensor will be a linear displacement current, a linear component of gating 
current and it would not be present in the reported gating current determined by 
some linear subtraction procedures. In particular, if the arginine system is 
present at the 'control' potential contributing a current linear in potential, this 
problem would occur. Of course, if the arginine system is immobilized and 
inactivated at the control potential and so contributes no current flow under that 
condition, this problem would not occur. 

Other systems may contribute to the linear displacement current as well, for 
example, i) all sorts of experimental and instrumentation artifacts and ii) 
displacement current in the conduction channel itself. The conduction channel of 
field effect transistors produces a large displacement current often characterized 
as a capacitance that involves diffusion and is described by drift diffusion 
equations quite similar to the PNP equations of the open conduction channel. 

Most systems have substantial motions that are linear in voltage (even if the 
system is labeled 'nonlinear'). The linear term is present in most systems, just as 
it is present in most Taylor expansions of nonlinear functions. 

The linear component that can be missed in experiments, and removed in 
these calculations, may have functional and structural significance. The voltage 
sensor works by sensing voltage, for example, by producing a motion of arginines. 
That motion—the response of the voltage sensor in this model—includes a linear 
component. The signal passed to the conduction channel, to control gating, is 
likely to include or depend on the linear component of sensor function. Confusion 
will result if a significant linear component exists and is ignored when a model is 
created that links the voltage sensor to the gating process of the conduction 
channel. Direct measurements of the movement of arginines (e.g., with optical 
methods) are likely to include the linear component and so should not agree with 
experimental measurements of gating current or with the currents reported here 
if the linear component exists and is significant in size. 

If the P/4 procedure subtracts a charge movement in a control system in 
which the arginines do not move at all (because they are immobilized and 
inactivated, in that sense), then the resulting estimate of gating current will 
contain a component linear in voltage. Thus, the interpretation of the corrected 
record depends on the details of immobilization and inactivation, topics that are 
beyond the scope of this paper and our present work. 

REFERENCES 

[1] Horng, T.L., T.C. Lin, C. Liu, and B. Eisenberg. 2012. PNP equations with steric 


9 



effects: a model of ion flow through channels. /. Phys. Chem. B. 
116 (3 7): 1142 2-11441. 

[2] Lin, T.C., and B. Eisenberg. 2014. A new approach to the Lennard-Jones 
potential and a new model: PNP-steric equations. Comm. Math. Sci. 
12(1):149-173. 

[3] Eisenberg, B., Y.K. Hyon, and C. Liu. 2010. Energy variational analysis EnVarA 
of ions in water and channels: field theory for primitive models of complex ionic 
fluids./. Chem. Phys. 133:104104. 

[4] Bezanilla, F., E. Perozo, and E. Stefani. 1994. Gating of K + channels : II. The 
components of gating currents and a model of channel. Biophy. J. 66:1011-1021. 

[5] Eisenberg, B. 2016a. Conservation of Current and Conservation of Charge. 
https://arxiv.org/abs/1609.09175. 

[6] Eisenberg, B. 2016b. Maxwell Matters, https://arxiv.org/pdf/1607.06691 . 

[7] Eisenberg, B., X. Oriols, and D. Ferry. 2017. Dynamics of current, charge, and 
mass. Mol. Based Math. Biol. 5:78-115. https://arxiv.org/abs/1708.0740Q . 

[8] Trefethen, L.N. 2000. Spectral Methods in MATLAB. SIAM, Philadelphia. 

[9] Ascher, U.M., and L.R. Petzold. 1998. Computer methods for ordinary 
differential equations and differential-algebraic equations. SIAM, Philadelphia. 

[10] Shampine, L.F., and M.W. Reichelt. 1997. The MATLAB ODE Suite. SIAM J. Sci. 
Comput. 18:1-22. 

[11] Shampine, L.F., M.W. Reichelt, and ].A. Kierzenka. 1999. Solving Index-1 DAEs 
in MATLAB and Simulink. SIAM Rev. 41:538-552. 

[12] Sigg, D., F. Bezanilla, and E. Stefani. 2003. Fast gating in the Shaker K+ 
channel and the energy landscape of activation. Proc. Natl. Acad. Sci. USA. 
100:7611-7615. 

[13] Stefani, E., and F. Bezanilla. 1997. Voltage dependence of the early events in 
voltage gating. Biophys.J. 72:131. 

[14] Stefani, E., D. Sigg, and F. Bezanilla. 2000. Correlation between the early 
component of gating current and total gating current in Shaker K channels. 
Biophys.J. 78:7. 

[15] Forster, I.C., and N.G. Greeff. 1992. The early phase of sodium channel gating 
current in the squid giant axon. Characteristics of a fast component of 
displacement charge movement. Eur. Biophys.J. 21(2): 99-116. 

[16] Schneider, M.F., and W.K. Chandler. 1973. Voltage dependent charge 
movement in skeletal muscle: a possible step in excitation contraction coupling. 
Nature. 242:244-246. 


10 




[17] Armstrong, C.M., and F. Bezanilla. 1974. Charge movement associated with 
the opening and closing of the activation gates of the sodium channels. /. Gen. 
Physiol 63:533-552. 

[18] Bezanilla, F., and J. Vergara. 1980. Properties of Excitable membranes. In 
Membrane structure and Function. E.E. Bittar, editor. Vol. II, Ch. 2. J. Wiley & Sons, 
New York, 53-113. 

[19] Bezanilla, F., and C.M. Armstrong. 1977. Inactivation of the sodium channel. I. 
Sodium current experiments./. Gen. Physiol. 70(5):549-566. 

[20] Fernandez, J.M., F. Bezanilla, and R.E. Taylor. 1982. Distribution and kinetics 
of membrane dielectric polarization. II. Frequency domain studies of gating 
currents./. Gen. Physiol. 79(l):41-67. 


11 



MATLAB CODE 

function arginine_s4qsbias 
global zl z2 z3 N1 N2 N3 Dzl Dz2 Dz3 
global DAzl DAz2 DAz3 Azl Az2 Az3 
global sphi K KS b S4offset LR L Qint 
global NaL NaR C1L C1R phiL phiR phiLlow 
global DNa Dq Dq2 DC1 

global zNa zCl zq Lap Gamma 1 Gamma2 Gamma3 
global iplot 

global ton toff Vwidtht 
global gij Wz 
format short e 

iclose=input('Input 1 to close old figures, default is 1: '); 
if isempty(iclose) 
iclose=l; 
end 

if iclose==l 
close all; 
end 

icon=input('Input 1 for preparing initial condition, 2 for full simulation, else for seeing old results, default is 1: '); 

if isempty(icon) 

icon=l; 

end 

if icon==l I icon==2 
if icon==2 

dir arginine-s4qsbias*.mat 

dataini=input('Input the file to load for initial condition: 
dataini=dataini(dataini~=''); 
tit2=['load' dataini]; 
eval(tit2); 

uO=uall(end,:); uO=uO(:); 
phiLlow=phiL; 

phiL=input('Input command voltage phiL, default is 0: '); 
if isempty(phiL) 
phiL=0; 
end 

ton=input('Input on time of phiL, default is 10: '); 
if isempty(ton) 
ton=10; 
end 

toff=input('Input off time of phiL, default is 150: '); 
if isempty(toff) 
toff=150; 
end 

Vwidtht=input('Input rise/fall rate of phiL, default is 5: '); 
if isempty (Vwidtht) 

Vwidtht=5; 

end 

Dq 

Dq=input('Input diffusion coefficient of arginine inside antechamber, default is 50: '); 
if isempty(Dq) 

Dq=50; 

end 

Dq2 
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Dq2=input('Input diffusion coefficient of arginine inside channel, default is 50: '); 
if isempty(Dq2) 

Dq2=50; 

end 

tspan=input('Input time span of simulation, default is 0:0.01:300: '); 

if isempty(tspan) 

tspan=(0:0.01:300)'; 

end 

else 

% valence: 

% zNa=input('Input zNa, default is 1: '); 

% if isempty(zNa) 
zNa=l; 

% end 

% zCl=input('Input zCl, default is -1:'); 

% if isempty(zCl) 
zCl=-l; 

% end 
% grid: 

Nzl=input('Input number of intervals in zone 1, default is 70: '); 
if isempty(Nzl) 

Nzl=70; 

end 

Nl=Nzl+l; Nz3=Nzl; N3=N1; 

Nz2=input('Input number of intervals in zone 2, default is 50: '); 
if isempty(Nz2) 

Nz2=50; 

end 

N2=Nz2+l; 

number_of_grids=[N 1 N2 N3] 

% geometry: 

L=0.7; % channel length with characteristic length lnm 

LR= 1.5; % antechamber length 

AL=0.15; % channel radius 

AR= 1; % antechamber radius 

Geometry_data=[L LR AL AR] 

% bulk diffusion coefficient: Na 1.33e-9; Cl 2.03e-9; arginine 0.7e-9 in m A 2/s 
DNa=l; DC1=1.53; % false diffusion coefficient for DNa and DC1 

Dq=input('Input diffusion coefficient of arginine inside antechamber, default is 5: '); 
if isempty(Dq) 

Dq=5; 

end 

Dq2=input('Input diffusion coefficient of arginine inside channel, default is 5: '); 
if isempty(Dq2) 

Dq2=5; 

end 

Diffusion_coefficient=[DNa DC1 Dq Dq2] 

% Gamma: 

Gamma l=input('Input Gamma in zone 1, default is 1: '); 
if isempty(Gammal) 

Gamma 1=1; 
end 

Gamma3=Gammal; 

Gamma2=input('Input Gamma in zone 2, default is 0.1:'); 
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if isempty(Gamma2) 

Gamma2=0.1; 

end 

Gamma_data=[Gammal Gamma2 Gamma3] 

% Boundary conditions: 

NaL= 1 ;NaR= 1 ;ClL=NaL;ClR=NaR; 

phiL=input('Input phiL to reach for its steady state, default is -3.6: '); 
if isempty(phiL) 
phiL=-3.6; 
end 

phiR=0; 

Boundary_condition_data=[NaL NaR C1L C1R phiL phiR] 

% collocation matrix 

% Left reservoir length LR, radius AR 

icho=2; 

if icho==l 

[Dxi,xi] =cheb(Nz 1); 

else 

[xi,tmp] = chebdif(Nzl+l,l); Dxi=reshape(tmp,Nzl+l,Nzl+l); 
end 

Dxi=fliplr(flipud(Dxi)); xi=flipud(xi); zl=(xi+l)*LR/2; Dzl=Dxi*2/LR; Dzzl=Dzl A 2; 
eyezl=eye(Nl); gzl=AR*ones(size(zl)); Azl=pi*gzl. A 2; DAzl=diag(l./Azl)*Dzl; 

% port length L, radius AL 
if icho==l 
[Dxi,xi] =cheb(Nz2); 
else 

[xi,tmp] = chebdif(Nz2+l,l); Dxi=reshape(tmp,Nz2+l,Nz2+l); 
end 

Dxi=fliplr(flipud(Dxi)); xi=flipud(xi); stp=LR; z2=(xi+l)*L/2+stp; Dz2=Dxi*2/L; Dzz2=Dz2 A 2; 
eyez2=eye(N2); gz2=AL*ones(size(z2)); Az2=pi*gz2. A 2; DAz2=diag(l./Az2)*Dz2; 

% right reservoir length LR, radius AR 

if icho==l 

[Dxi,xi] =cheb(Nz3); 

else 

[xi,tmp] = chebdif(Nz3+l,l); Dxi=reshape(tmp,Nz3+l,Nz3+l); 
end 

Dxi=fliplr(flipud(Dxi)); xi=flipud(xi); stp=LR+L; z3=(xi+l)*LR/2+stp; Dz3=Dxi*2/LR; Dzz3=Dz3 A 2; 
eyez3=eye(N3); gz3=AR*ones(size(z3)); Az3=pi*gz3. A 2; DAz3=diag(l./Az3)*Dz3; 

% figure(l); plot(zl,gzl,'b','LineWidth',2); ylim([-8 8]); grid on; xlabel('z','FontSize',18); ylabel('r','FontSize',18); hold 
on; 

% plot(z2,gz2,'b','LineWidth',2); plot(z3,gz3,'b','LineWidth',2); 

% plot(zl,-gzl,'b','LineWidth',2); plot(z2,-gz2,'b','LineWidth',2); plot(z3,-gz3,'b','LineWidth',2); 

% plot([LR LR],[-AR AR],'r','LineWidth',2); plot([LR+L LR+L],[-AR AR],'r','LineWidth',2); hold off; 

% title('3-zone channel shape','FontSize',18); drawnow; 

% Laplace operator for electric potential: 

Lapl=Dzzl; Lapl(l,:)=eyezl(l,:); Lapl(end,:)=Gamrnal*Azl(l)*Dzl(end,:); 

Lap2=Dzz2; Lap2(l,:)=eyez2(l,:); Lap2(end,:)=Gamma2*Az2(l)*Dz2(end,:); 

Lap3=Dzz3; Lap3(l,:)=eyez3(l,:); Lap3(end,:)=eyez3(end,:); 

Lap=blkdiag(Lapl,Lap2,Lap3); 

Lap(Nl,Nl+l:Nl-t-N2)=-Gamma2*Az2(l)*Dz2(l,:); Lap(Nl+l,l:Nl)=-eyezl(end,:); 
Lap(Nl+N2,Nl+N2+l:end)=-Gamma3*Az3(l)*Dz3(l,:); Lap(Nl+N2+l,Nl+l:Nl+N2)=-eyez2(end,:); 

% mass matrix: Na, Cl, 4 arginine: 

massl=ones(size(zl)); massl([l end])=0; masslnull=zeros(size(zl)); 
mass2=ones(size(z2)); mass2([l end])=0; 
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mass3=ones(size(z3)); mass3([l end])=0; mass3null=zeros(size(z3)); 
M=sparse(diag([masslnull(:);masslnull(:);massl(:);massl(:);massl(:);massl(:); ... 
mass2(:);mass2(:);mass2(:);mass2(:);... 

mass3null(:);mass3null(:);mass3(:);mass3(:);mass3(:);mass3(:);l])); 

% initial condition: 

Na01=NaL*ones(size(zl)); Na03=NaR*ones(size(z3)); C101=Na01; C103=Na03; 

% Arginine initial distribution 

Q=input('Input initial uniform concentration distribution for each arginine at left antechamber, default is 1/10: '); 
if isempty(Q) 

Q=l/10; 

end 

disp('Volume of each arginine calcualted from its initial concentration distribution:'); 

Qint=Az 1 (1) * LR * Q 
zq=input('Input zq, default is 1: '); 
if isempty(zq) 

zq=l; % valence of arginine 
end 

% if abs(phiL+1.2)>=2 

% filterl=(l-tanh(10*(zl-LR-0.2)))/2; filter2=(l-tanh(10*(z2-LR-0.2)))/2; filter3=(l-tanh(10*(z3-LR-0.2)))/2; 

% q01=Q*filterl; q02=Q*filter2; q03=Q*filter3; 

% else 

Quniform=Qint/( Az2( 1) *L+Az 1(1) *LR+Az3( 1) *LR); 

q01=ones(size(zl))*Quniform; q02=ones(size(z2))*Quniform; q03=ones(size(z3))*Quniform; 

% end 

% figure(lOO); plot(zl,q01,'b','LineWidth',2); xlabel('z','FontSize',18); ylabel('initial arginine','FontSize',18); grid on; 
hold on; 

% plot(z2,q02,'g','LineWidth',2); plot(z3,q03,'rVLineWidth',2); ylim([-0.1*Q 1.1*Q]); hold off; drawnow; 
u0= 

[Na01(:);C101(:);q01(:);q01(:);q01(:);q01(:);q02(:);q02(:);q02(:);q02(:);Na03(:);C103(:);q03(:);q03(:);q03(:);q03(:);LR+ 

L/2]; 

% steric effect: 

gij=input('Input steric effect parameter g, default is 0.5: '); 
if isempty(gij) 
gij=0.5; 
end 

% potential trap: 

sphi=0.2; % half of the interval between spring anchoring position on S4 for each arginine 

K=input('Input spring constant K for each arginine connected to S4, default is 3:'); 
if isempty(K) 

K=3; 

end 

KS=input('Input spring constant for S4, needing to be larger than 4K, default is 3:'); 
if isempty(KS) 

KS=3; 

end 

b=input('Input damping coefficient for S4, must be positive, too small will make equations stiff, default is 1.5:'); 
if isempty(b) 
b=1.5; 
end 

S4offset=input('Input offset of equilibrium position of S4 from middle of channel, must be positive, default is 1.591: '); 
if isempty(S4offset) 

S4offset= 1.591; 
end 

Spring_data=[sphi K KS b S4offset] 
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% energy barrier inside channel: 

Wmag=input('Input magnitude of energy barrier V inside channel, default is 5: '); 
if isempty(Wmag) 

Wmag=5; 

end 

Wwidth=input('Input rise/fall rate in space for energy barrier V above, the purpose is to smooth, default is 5: '); 
if isempty(Wwidth) 

Wwidth=5; 

end 

energy_barrier_data=[Wmag Wwidth] 

W=Wmag*tanh(Wwidth*(z2-LR))-Wmag*tanh(Wwidth*(z2-LR-L))-Wmag; 

Wz=Wmag*Wwidth*sech(Wwidth*(z2-LR)). A 2-Wmag*Wwidth*sech(Wwidth*(z2-LR-L)). A 2; 

% figure(2); subplot(l,2,l); plot(z2,W,'b','LineWidth',2); xlabel('z','FontSize',18); ylabel('energy barrier 
V','FontSize',18); grid on; 

% subplot(l,2,2); plot(z2,Wz,'b','LineWidth',2); xlabel('z','FontSize',18); ylabel('dV/dz','FontSize',18); grid on; 
drawnow; 

% ODE parameters: 

% if (abs(phiL+1.2)<=l) 

tspan=input('Input time span of calculation to reach steady state, default is 0:0.01:800:'); 

if isempty(tspan) 

tspan=(0:0.01:800)'; 

end 

% else 

% tspan=input('Input time span of calculation to reach steady state, default is 0:0.01:400: '); 

% if isempty(tspan) 

% tspan=(0:0.01:400)'; 

% end 
% end 
end 

mxstep=input('Input max numerical time step allowed, default is 0.005: '); 

if isempty(mxstep) 

mxstep=0.005; 

end 

iplot=input('Input 1 to plot for the purpose of debugging during computation, default is 0:'); 

if isempty(iplot) 

iplot=0; 

end 

options=odeset('RelTol', le-3,'AbsTol', le-5,'Mass',M,'MaxStep',mxstep); 

%tic; 

ton 

% pause; 

[t,uall] =ode 15 s( @ pnp 1 d,tspan,u0,options); 

%et=toc; 
if -isempty(ton) 

tit=['arginine-s4qsbias-gij- num2str(gij) '-phiLlow=' num2str(phiLlow) '-phiL=' num2str(phiL) '-K- num2str(K) '-KS 
num2str(KS)... 

'-b=' num2str(b) '-S4offset=' num2str(S4offset) '-Dargl=' num2str(Dq) '-Darg2=' num2str(Dq2) '-Q=' num2str(Q) '- 
Gamma2=' num2str(Gamma2)... 

'-Wmag=' num2str(Wmag) '-maxstep=' num2str(mxstep) '-tend- num2str(t(end)) '-zq=' num2str(zq) '-Vwidtht=' 
num2str(Vwidtht)] 

% tit=['arginine-s4qsbias-gij- num2str(gij) '-phiLlow- num2str(phiLlow) '-phiL=' num2str(phiL) '-K— num2str(K) '- 
KS=' num2str(KS) ... 

% '-b=' num2str(b) '-S4offset=' num2str(S4offset) '-Dargl=' num2str(Dq) '-Darg2=' num2str(Dq2) '-Q=' num2str(Q) 

Gamma2=' num2str(Gamma2)... 
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% '-Wmag=' num2str(Wmag) '-maxstep=' num2str(mxstep) '-toff=' num2str(toff)] 
else 

tit=['arginine-s4qsbias-gij=' num2str(gij) '-phiL=' num2str(phiL) '-K=' num2str(K) '-KS=' num2str(KS)... 

'-b=' num2str(b) '-S4offset=' num2str(S4offset) '-Dargl=' num2str(Dq) '-Darg2=' num2str(Dq2) '-Q=' num2str(Q) '- 
Gamma2=' num2str(Gamma2)... 

'-Wmag=' num2str(Wmag) '-maxstep=' num2str(mxstep) '-tend- num2str(t(end)) '-zq=' num2str(zq)] 
end 

titl=['save' tit '.mat']; 

eval(titl); 

else 

close all; 

dir arginine-s4qsbias*.mat 

dataini=input('Input the file to load: ','s'); 

dataini=dataini(dataini~=''); 

tit2=['load' dataini]; 

eval(tit2); 

end 

Il=[]; I3=[]; 

I2a=[];I2b=[];I2c=[];I2d=[]; 

Ql=[];Q2=[];Q3=[];Q2half=[]; 
phiLv=[]; tv=[]; 

Nalall=[]; Cllall=[]; Na3all=[]; C13all=[]; NamCll=[]; NamC13=[]; 

qlaall=[]; qlball=[]; qlcall=[]; qldall=[]; 

q2aall=[]; q2ball=[]; q2call=[]; q2dall=[]; 

q3aall=[]; q3ball=[]; q3call=[]; q3dall=[]; 

philall=[]; phi2all=[]; phi3all=[]; 

S4dispall=[]; qaposall=[]; qbposall=[]; qcposall=[]; qdposall=[]; 

int=input('Input output time interval, default is 20: '); 

if isempty(int) 

int=20; 

end 

for i=l:int:length(t) 
tv=[tv;t(i)]; 
u=uall(i,:); u=u.'; 

Nal=u(l:Nl); Cll=u(Nl+l:2*Nl); qla=u(2*Nl+l:3*Nl); qlb=u(3*Nl+l:4*Nl); qlc=u(4*Nl+l:5*Nl); 
qld=u(5*Nl+l:6*Nl); stp=6*Nl; 

q2a=u(stp+l:stp+N2); q2b=u(stp+N2+l:stp+2*N2); q2c=u(stp+2*N2+l:stp+3*N2); q2d=u(stp+3*N2+l:stp+4*N2); 
stp=6*Nl+4*N2; 

Na3=u(stp+l:stp+N3); C13=u(stp+N3+l:stp+2*N3); q3a=u(stp+2*N3+l:stp+3*N3); q3b=u(stp+3*N3+l:stp+4*N3); 
q3c=u(stp+4*N3+l:stp+5*N3); q3d=u(stp+5*N3+l:stp+6 :|: N3); S4disp=u(end); 

Nalall=[Nalall;Nal.']; Cl 1 all=[Cl 1 all;Cl 1.']; Na3all=[Na3all;Na3.']; C13all=[C13all;C13.']; 
qlaall=[qlaall;qla.']; qlball=[qlball;qlb.']; qlcall=[qlcall;qlc.']; qldall=[qldall;qld.']; 
q2aall=[q2aall;q2a.']; q2ball=[q2ball;q2b.']; q2call=[q2call;q2c.']; q2dall=[q2dall;q2d.']; 
q3aall=[q3aall;q3a.']; q3ball=[q3ball;q3b.']; q3call=[q3call;q3c.']; q3dall=[q3dall;q3d.']; 

ql=qla+qlb+qlc+qld; warning off; tmp=Dzl\(Azl.*ql); warning on; qlint=tmp(end)-tmp(l); Ql=[Ql;qlint]; 
q2=q2a+q2b+q2c+q2d; warning off; tmp=Dz2\(Az2.*q2); warning on; q2int=tmp(end)-tmp(l); Q2=[Q2;q2int]; 
q2half=q2a+q2b+q2c+q2d; q2half(z2>LR+L/2)=0; warning off; tmp=Dz2\(Az2.*q2half); 
warning on; q2halfint=tmp(end)-tmp(l); Q2half=[Q2half;q2halfint]; 

q3=q3a+q3b+q3c+q3d; warning off; tmp=Dz3\(Az3.*q3); warning on; q3int=tmp(end)-tmp(l); Q3=[Q3;q3int]; 

warning off; tmp=Dzl\(Azl.*(Nal-Cll)); warning on; NamCllint=tmp(end)-tmp(l); NamCll=[NamCll;NamCllint]; 

warning off; tmp=Dz3\(Az3.*(Na3-C13)); warning on; NamC13int=tmp(end)-tmp(l); NamC13=[NamC13;NamC13int]; 

Vwidth=Vwidtht; 

if isempty(ton) 

phiLactual=phiL; 
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else 

phiLactual=((phiL-phiLlow)*tanh(Vwidth*(t(i)-(ton+l)))-(phiL-phiLlow)*tanh(Vwidth*(t(i)-(toff-l))))/2+phiLlow; 

end 

phiLv=[phiLv;phiLactual]; 

rhsl=-(zNa*Nal+zCl*Cll+zq*(qla+qlb+qlc+qld))/Gammal; rhsl(l)=phiLactual; rhsl(end)=0; 
rhs2=-(zq*(q2a+q2b+q2c+q2d))/Gamma2; rhs2([l end])=0; 

rhs3=-(zNa*Na3+zCl*C13+zq*(q3a+q3b+q3c+q3d))/Gamma3; rhs3(l)=0; rhs3(end)=phiR; 
rhs=[rhsl;rhs2;rhs3]; warning off; phi=Lap\rhs; warning on; 
phil=phi(l:Nl); phi2=phi(Nl+l:Nl+N2); phi3=phi(Nl+N2+l:end); 
phi 1 all=[phi 1 all;phi 1.']; phi2all=[phi2all;phi2.']; phi3all=[phi3all;phi3.']; 

Nazl=Dzl*Nal; Naz3=Dz3*Na3; Clzl=Dzl*Cll; Clz3=Dz3*C13; 
qzla=Dzl*qla; qzlb=Dzl*qlb; qzlc=Dzl*qlc; qzld=Dzl*qld; 
qz2a=Dz2*q2a; qz2b=Dz2*q2b; qz2c=Dz2*q2c; qz2d=Dz2*q2d; 
qz3a=Dz3*q3a; qz3b=Dz3*q3b; qz3c=Dz3*q3c; qz3d=Dz3*q3d; 
phizl=Dzl*phil; phiz3=Dz3*phi3; phiz2=Dz2*phi2; 
warning off; 

tmp=Azl(l)*qla.*zl; tmp=Dzl\tmp; tmpl=tmp(end)-tmp(l); 
tmp=Az2(l)*q2a.*z2; tmp=Dz2\tmp; tmp2=tmp(end)-tmp(l); 
tmp=Az3(l)*q3a.*z3; tmp=Dz3\tmp; tmp3=tmp(end)-tmp(l); 
qapos=(tmp 1 +tmp2+tmp3)/Qint; qaposall=[qaposall;qapos]; 
tmp=Azl(l)*qlb.*zl; tmp=Dzl\tmp; tmpl=tmp(end)-tmp(l); 
tmp=Az2(l)*q2b.*z2; tmp=Dz2\tmp; tmp2=tmp(end)-tmp(l); 
tmp=Az3(l)*q3b.*z3; tmp=Dz3\tmp; tmp3=tmp(end)-tmp(l); 
qbpos=(tmpl+tmp2+tmp3)/Qint; qbposall=[qbposall;qbpos]; 
tmp=Azl(l)*qlc.*zl; tmp=Dzl\tmp; tmpl=tmp(end)-tmp(l); 
tmp=Az2(l)*q2c.*z2; tmp=Dz2\tmp; tmp2=tmp(end)-tmp(l); 
tmp=Az3(l)*q3c.*z3; tmp=Dz3\tmp; tmp3=tmp(end)-tmp(l); 
qcpos=(tmp 1 +tmp2+tmp3)/Qint; qcposall=[qcposall;qcpos]; 
tmp=Azl(l)*qld.*zl; tmp=Dzl\tmp; tmpl=tmp(end)-tmp(l); 
tmp=Az2(l)*q2d.*z2; tmp=Dz2\tmp; tmp2=tmp(end)-tmp(l); 
tmp=Az3(l)*q3d.*z3; tmp=Dz3\tmp; tmp3=tmp(end)-tmp(l); 
qdpos=(tmpl+tmp2+tmp3)/Qint; qdposall=[qdposall;qdpos]; 
warning on; 

S4dispall=[S4dispall;S4disp]; 

Vzla=K*(zl-(S4disp-3*sphi)); Vzlb=K*(zl-(S4disp-sphi)); Vzlc=K*(zl-(S4disp+sphi)); Vzld=K*(zl- 
(S4disp+3*sphi)); 

Vz2a=K*(z2-(S4disp-3*sphi)); Vz2b=K*(z2-(S4disp-sphi)); Vz2c=K*(z2-(S4disp+sphi)); Vz2d=K*(z2- 
(S4disp+3*sphi)); 

Vz3a=K*(z3-(S4disp-3*sphi)); Vz3b=K*(z3-(S4disp-sphi)); Vz3c=K*(z3-(S4disp+sphi)); Vz3d=K*(z3- 
(S4disp+3*sphi)); 

JNazl=-Azl.*DNa.*(Nazl+zNa*Nal.*phizl); 

JClz 1=-Az 1. *DC1. *(Clz 1 +zCl*Cl 1. *phiz 1); 

Jqzla=-Azl.*Dq.*(qzla+zq*qla.*phizl+qla.*Vzla+gij*qla.*(qzlb+qzlc+qzld)); 

Jqzlb=-Azl.*Dq.*(qzlb+zq*qlb.*phizl+qlb.*Vzlb+gij*qlb.*(qzla+qzlc+qzld)); 

Jqzlc=-Azl.*Dq.*(qzlc+zq*qlc.*phizl+qlc.*Vzlc+gij*qlc.*(qzla+qzlb+qzld)); 

Jqz 1 d=-Az 1. *Dq. *(qz 1 d+zq*q 1 d. *phiz l+qld.*Vzl d+gij *q 1 d. *(qz 1 a+qz 1 b+qz 1 c)); 

Jqz2a=-Az2.*Dq2.*(qz2a+zq*q2a.*phiz2+q2a.*(Wz+Vz2a)+gij*q2a.*(qz2b+qz2c+qz2d)); 

Jqz2b=-Az2.*Dq2.*(qz2b+zq*q2b.*phiz2+q2b.*(Wz+Vz2b)+gij*q2b.*(qz2a+qz2c+qz2d)); 

Jqz2c=-Az2.*Dq2.*(qz2c+zq*q2c.*phiz2+q2c.*(Wz+Vz2c)+gij*q2c.*(qz2a+qz2b+qz2d)); 

Jqz2d=-Az2.*Dq2.*(qz2d+zq*q2d.*phiz2+q2d.*(Wz+Vz2d)+gij*q2d.*(qz2a+qz2b+qz2c)); 

JNaz3=-Az3.*DNa.*(Naz3+zNa*Na3.*phiz3); 

JClz3=-Az3.*DCl.*(Clz3+zCl*C13.*phiz3); 

Jqz3a=-Az3.*Dq.*(qz3a+zq*q3a.*phiz3+q3a.*Vz3a+gij*q3a.*(qz3b+qz3c+qz3d)); 

Jqz3b=-Az3.*Dq.*(qz3b+zq*q3b.*phiz3+q3b.*Vz3b+gij*q3b.*(qz3a+qz3c+qz3d)); 
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Jqz3c=-Az3.*Dq.*(qz3c+zq*q3c.*phiz3+q3c.*Vz3c+gij*q3c.*(qz3a+qz3b+qz3d)); 

Jqz3d=-Az3.*Dq.*(qz3d+zq*q3d.*phiz3+q3d.*Vz3d+gij*q3d.*(qz3a+qz3b+qz3c)); 

Jl=zNa*JNazl+zCl*JClzl+zq*(Jqzla+Jqzlb+Jqzlc+Jqzld); 

J3=zNa*JNaz3+zCl*JClz3+zq*(Jqz3a+Jqz3b+Jqz3c+Jqz3d); 

J2a=zq*Jqz2a; J2b=zq*Jqz2b; J2c=zq*Jqz2c; J2d=zq*Jqz2d; 

I1=[I1 ;J 1.']; I3=[I3;J3.']; 

I2a=[I2a;J2a.']; I2b=[I2b;J2b.']; I2c=[I2c;J2c.']; I2d=[I2d;J2d.']; 
end 

I2=I2a+I2b+I2c+I2d; 

dt=tv(2)-tv(l); 

% figure(4); 

% subplot(l,2,l); plot(zl,gzl,'b','LineWidth',2); ylim([-8 8]); grid on; xlabel('z','FontSize',18); ylabel('r','FontSize',18); 
hold on; 

% plot(z2,gz2,'b','LineWidth',2); plot(z3,gz3,'b','LineWidth',2); 

% plot(zl,-gzl,'b','LineWidth',2); plot(z2,-gz2,'b','LineWidth',2); plot(z3,-gz3,'b','LineWidth',2); 

% plot([LR LR],[-AR AR],'r','LineWidth',2); plot([LR+L LR+L],[-AR AR],'r','LineWidth',2); hold off; 

% title('3-zone channel shape','FontSize',18); drawnow; 

% subplot(l,2,2); plot(z2,W,'b','LineWidth',2); xlabel('z','FontSize',18); ylabel('energy barrier V','FontSize',18); 

% grid on; xlim([LR L+LR]); ylim([-l Wmag+1]); 

% I at several location 
[tmp,izmid]=min(abs(z2-(LR+L/2))); 

Imid=I2(:,izmid); Imida=I2a(:,izmid); Imidb=I2b(:,izmid); Imidc=I2c(:,izmid); Imidd=I2d(:,izmid); 

Imidtmp=Imid( 1 :end-1); 

NImid=length(Imidtmp); Ttv=tv(end); 

ky=[0:NImid/2 -NImid/2+l:-l]'; iky=li*[eps:NImid/2-l eps -NImid/2+l:-l]'; 

Imidh=fft(Imidtmp); tmp=Imidh./iky; tmp(l)=0; tmp(NImid/2+l)=0; Imidint=real(ifft(tmp)*Ttv/(2*pi)); Imidint= 
[Imidint;Imidint(end)]; 

Imidint=Imidint-Imidint( 1); 

Qmid=zq*(Ql+Q2half)+NamCll; ImidQ=zeros(size(tv)); ImidQ(2:end-l)=(Qmid(3:end)-Qmid(l:end-2))/(2*dt); 
ImidQ(l)=(-3*Qmid(l)+4*Qmid(2)-Qmid(3))/(2*dt); ImidQ(end)=(-3*Qmid(end)+4*Qmid(end-l)-Qmid(end- 
2))/(-2*dt); 

ImidQ=(Az2( 1 )*Imid-ImidQ)/Az 1(1); 
if -isempty(ton) 
figure (7); 

subplot(3,l,l); plot(tv,phiLv*25,'ColorVb','LineWidth',2); xlabel('t (a.u.)','FontSize',20); 
ylabel('V (mV)','FontSize',20); grid on; xlim([0 tv(end)]); 
ylim( [(phiLlow-1) *25 (phiL+1) *25 ]); 
subplot(3,l,2); 

plot(tv,Imid,'Color','b','LineWidth',2); xlabel('t (a.u.)','FontSize',20); ylabel('I (a.u.)','FontSize',20); grid on; xlim([0 
tv(end)]); 

hold on; plot(tv,Imida,'Color','r','LineWidth',2); plot(tv,Imidb,'Color','k','LineWidth',2); 
plot(tv,Imidc,'Color','m','LineWidth',2); plot(tv,Imidd,'Color','c','LineWidth',2); hold off; 
hh=legend('total',T,'2','3','4'); 
subplot(3,l,3); 

plot(tv,ImidQ,'Color','b','LineWidth',2); xlabel('t (a.u.)','FontSize',20); ylabel('I voltage clamp (a.u.)','FontSize',20); 

grid on; xlim([0 tv(end)]); 

end 

if -isempty(ton) 

[tmp ,imax] =max(Imid); 

[tmp,imin]=min(abs(tv-140)); 
xdata=tv(imax: imin); ydata=Imid(imax: imin); 
options=optimset('MaxFunEvals',10000,'MaxIter', 10000); 
xO = [0.02;10;0.02;10]; 

Ib=[0;0;0;0];ub=[le4;le4;le4;le4]; 
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[x2,resnorm2] = Isqcurvefit(@myfun2,x0,xdata,ydatajb,ub,options); 
x2 

resnorm2 

ydatafit2=myfun2(x2,xdata); 
xO = [0.02;10;0.02;10;0.02;10]; 
lb=[0;0;0;0;0;0] ;ub=[ le4; le4; le4; le4; le4; le4]; 

[x3,resnorm3] = Isqcurvefit(@myfun3,x0,xdata,ydatajb,ub,options); 
x3 

resnorm3 

ydatafit3=myfun3(x3 ,xdata); 

xO = [0.02;10;0.02;10;0.02;10;0.02;10]; 

lb=[0;0;0;0;0;0;0;0] ;ub=[ le4; le4; le4; le4; le4; le4; le4; le4]; 

[x4,resnorm4] = Isqcurvefit(@myfun4,x0,xdata,ydata,lb,ub,options); 
x4 

resnorm4 

ydatafit4=myfun4(x4,xdata); 
figure (70); 

plot(xdata,ydata,'Color','b','lineWidth',2); xlabel('f,'FontSize',18); ylabel('Imiddle','FontSize',18); grid on; hold on; 
plot(xdata,ydatafit2,'ColorVg','lineWidth',2); plot(xdata,ydatafit3,'ColorVr','lineWidth',2); 
plot(xdata,ydatafit4,'ColorVkVlineWidth',2); hold off; legend('originar,'2','3','4'); 
xO = [0.02;10;0.02;10;1;1]; 
lb=[0;0;0;0;0;0] ;ub=[ le4; le4; le4; le4; le4; le4]; 

[x5,resnorm5] = Isqcurvefit(@myfun5,x0,xdata,ydatajb,ub,options); 
x5 

resnorm5 

ydatafit5 =myfun5 (x5 ,xdata); 

xO = [0.02;10;0.02;10;0.02;10;1;1;1]; 

lb=[0;0;0;0;0;0;0;0;0] ;ub=[ le4; le4; le4; le4; le4; le4; le4; le4; le4]; 

[x6,resnorm6] = Isqcurvefit(@myfun6,x0,xdata,ydatajb,ub,options); 
x6 

resnorm6 

ydatafit6=myfun6(x6,xdata); 

xO = [0.02;10;0.02;10;0.02;10;0.02;10;1;1;1;1]; 

Ib=[0;0;0;0;0;0;0;0;0;0;0;0];ub=[le4;le4;le4;le4;le4;le4;le4;le4;le4;le4;le4;le4]; 

[x7,resnorm7] = Isqcurvefit(@myfun7,x0,xdata,ydatajb,ub, options); 
x7 

resnorm7 

ydatafit7=myfun7 (x7 ,xdata); 
figure(71); 

plot(xdata,ydata,'Color','bVlineWidth',2); xlabel('f,'FontSize',18); ylabel('Imiddle','FontSize',18); grid on; hold on; 
plot(xdata,ydatafit5,'Color','g','lineWidth',2); plot(xdata,ydatafit6,'ColorVr','lineWidth',2); 
plot(xdata,ydatafit7,'ColorVkVlineWidth',2); hold off; legend('originar,'2','3','4'); 

Qint 

end 

V=Qint*4; 

figure(8); 

plot(tv,4*Ql/V,'Color','bVLineWidth',2); xlabel('t (a.u.)','FontSize\20); ylabel('arginine','FontSize',20); grid on; hold on 

plot(tv,4*Q3/V,'ColorVg','LineWidth',2); plot(tv,4*Q2/V,'ColorVr','LineWidth',2); 

%plot(tv,(Ql+Q2+Q3)/V,'Color';k','LineWidth',2); 

hold off; 

ylim([-0.5 4.5]); 

if isempty(ton) 

title(['Final Ql/V=' num2str(Ql(end)/V)Q3/V=' num2str(Q3(end)/V)]); 
end 
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hh=legend('zone l','zone 3','zone 2'); 
figure(81); 

plot(tv,zq*Ql+NamCll,'Color','bVLineWidth',2); xlabel('t (a.u.)','FontSize',20); ylabel('total net charge','FontSize',20); 
grid on; hold on; plot(tv,zq*Q3+NamC13,'Color','g','LineWidth',2); hold off; 
hh=legend('zone l','zone 3'); 
figure(9); 

% qpos=0.25*(qaposall+qbposall+qcposall+qdposall); 

% qposmax=max( [qaposall-qaposall( 1) ;qbposall-qbposall( 1) ;qcposall-qcposall( 1) ;qdposall-qdposall( 1)]); 

% Imidintmax=max(Imidint);Imidint=Imidint*qposmax/Imidintmax; 
subplot(2,l,l); 

plot(tvJmidint*4/V,'bVLineWidth',2); grid on; xlabel('t (a.u.)','FontSize',20); ylabel('charges moved, Q','FontSize',20); 
subplot(2,l,2); 

plot(tv,qaposall-qaposall(l),'rVLineWidth',2); grid on; xlabel('t (a.u.)','FontSize',20); 

ylabel('\Delta z_{i,CM}, \Delta z_{S4} (nm)','FontSize',20); hold on; 

plot(tv,qbposall-qbposall(l),'k','LineWidth',2); 

plot(tv,qcposall-qcposall(l),'m','LineWidth',2); 

plot(tv,qdposall-qdposall(l),'c','LineWidth',2); 

plot(tv,S4dispall-S4dispall(l),'b','LineWidth',2); 

hold off; 

hh=legend('\Delta z_{ l,CM}','\Delta z_{2,CM}','\Delta z_{3,CM}','\Delta z_{4,CM}','\Delta z_{S4}'); 

figure (91); 

subplot(2,2,l); 

% qpos=0.25*(qaposall+qbposall+qcposall+qdposall); 

% qposmax=max( [qaposall-qaposall( 1) ;qbposall-qbposall( 1) ;qcposall-qcposall( 1) ;qdposall-qdposall( 1)]); 

% Imidintmax=max(Imidint);Imidint=Imidint*qposmax/Imidintmax; 
plot(tv,qaposall,'r','LineWidth',2); grid on; xlabel('t (a.u.)','FontSize',20); 
ylabel('z_{i,CM}, z_{ S4} (nm)','FontSize',20); hold on; 
plot(tv,qbposall,'k','LineWidth',2); 
plot(tv,qcposall,'m','LineWidth',2); 
plot(tv,qdposall,'c','LineWidth',2); 
plot(tv,S4dispall,'b','LineWidth',2); 
hold off; 

hh=legend('z_ {1 ,CM} ','z_ {2,CM} ','z_ {3 ,CM} 7z_ {4,CM} ','z_{ S4}'); 
ymaxl=max([Nalall(:);Na3all(:);Cllall(:);C13all(:)]); 
y max2=max( [phi 1 all(:) ;phi2all(:) ;phi3 all(:) ]); 
ymin2=min( [phi 1 all(:) ;phi2all(:) ;phi3all(:)]); 

% ymax3=max(Ileft); 

% ymin3=min(Ileft); 
subplot(2,2,2); 

plot(tv,qaposall-qaposall(l),'rVLineWidth',2); grid on; xlabel('t (a.u.)','FontSize',20); 

ylabel('\Delta z_[i,CM}, \Delta z_[S4} (nm)','FontSize',20); hold on; 

plot(tv,qbposall-qbposall(l),'k','LineWidth',2); 

plot(tv,qcposall-qcposall(l),'m','LineWidth',2); 

plot(tv,qdposall-qdposall(l),'c','LineWidth',2); 

plot(tv,S4dispall-S4dispall(l),'b','LineWidth',2); 

hold off; 

subplot(2,2,3) 

qav=zeros(size(qaposall)); qbv=qav; qcv=qav; qdv=qav; S4v=qav; 

qav(2:end-l)=(qaposall(3:end)-qaposall(l:end-2))/(2*dt); 

qbv(2:end-l)=(qbposall(3:end)-qbposall(l:end-2))/(2*dt); 

qcv(2:end-l)=(qcposall(3:end)-qcposall(l:end-2))/(2*dt); 

qdv(2:end-l)=(qdposall(3:end)-qdposall(l:end-2))/(2*dt); 

S4v(2:end-l)=(S4dispall(3:end)-S4dispall(l:end-2))/(2*dt); 
qav( 1 )=qav(2); qav(end)=qav(end-1); 
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qbv( 1 )=qb v(2); qbv(end)=qbv(end-1); 
qcv( 1 )=qcv(2); qcv(end)=qcv(end-1); 
qdv( 1 )=qdv(2); qdv(end)=qdv(end-1); 

S4v( 1 )=S4v(2); S4v(end)=S4v(end-1); 

plot(tv,qav,'r','LineWidth',2); grid on; xlabel('t (a.u.)','FontSize',20); 

ylabel('v_{i,CM}, v_{S4} (nm)7FontSize',20); hold on; 

plot(tv,qbv,'k','LineWidth',2); 

plot(tv,qcv,'m','LineWidth',2); 

plot(tv ,qdv, 'c' ,'Line W idth',2); 

plot(tv,S4v,'b','LineWidth',2); 

subplot(2,2,4) 

plot(tv,qav,'r7LineWidth',2); grid on; xlabel('t (a.u.)','FontSize',20); 

ylabel('v_{i,CM}, v_{S4} (nm)','FontSize',20); hold on; 

plot(tv,qbv,'k','LineWidth',2); 

plot(tv,qcv,'m7LineWidth',2); 

plot(tv,qdv,'c','LineWidth',2); 

plot(tv,S4v,'b','LineWidth',2); 

if isempty(ton) 

figure(lO); 

plot(zl,Nalall(end,:)3,'b','LineWidth',2); hold on; 
plot(z 1 ,C11 all(end,:).','g','LineWidth',2); 
plot(z 1 ,q 1 aall(end,:).' ,'r', 'Line Width', 2); 
plot(zl,qlball(end,:).','k','LineWidth',2); 
plot(zl,qlc all(end,:). ’' ,'m', 'Line W idth', 2); 
plot(z 1 ,q ldall(end, :).','c','LineWidth',2); 
plot(z2 ,q2aall(end,:).' ,'r', 'Line W idth', 2); 
plot(z2,q2ball(end,:).','k','LineWidth',2); 
plot(z2,q2call(end,:).’,'m','LineWidth',2); 
plot(z2,q2dall(end,:).','c','LineWidth',2); 
plot(z3,q3aall(end,:).','r','LineWidth',2); 
plot(z3,q3ball(end,:).','k','LineWidth',2); 
plot(z3,q3call(end,:).’,'m','LineWidth',2); 
plot(z3,q3dall(end,:).','c','LineWidth',2); 
plot(z3,Na3all(end,:).','b','LineWidth',2); 
plot(z3,C13all(end,:).','g','Line Width',2); 

hold off; grid on; xlabel('z','FontSize',18); ylabel('Na, Cl, Arginine','FontSize',18); xlim([0 L+2*LR]); 

ylim([-0.1 ymaxl]); hh=legend('Na','Cl','c_l','c_2','c_3','c_4','Location','SouthEast'); 

else 

figure(lO); 

subplot(2,4,l); 

plot(zl,Nalall(l,:).','b','LineWidth',2); hold on; 

plot(z 1 ,C11 all( 1, :).','g','LineWidth',2); 

plot(zl,qlaall(l,:).' ,'r' ,'Line W idth', 2); 

plot(z 1 ,q 1 ball( 1,:).' ,'k', 'Line W idth', 2); 

plot(zl,qlcall(l,:).','m','LineWidth',2); 

plot(zl,qldall(l,:).','c','LineWidth',2); 

plot(z2,q2aall(l,:).','r','LineWidth',2); 

plot(z2,q2ball(l,:).','k','LineWidth',2); 

plot(z2,q2call(l,:).','m','LineWidth',2); 

plot(z2,q2dall(l,:).','c','LineWidth',2); 

plot(z3,q3aall(l,:).' ,'r' ,'Line W idth', 2); 

plot(z3 ,q3ball( 1,:).' ,'k', 'Line W idth', 2); 

plot(z3,q3call(l,:).','m','LineWidth',2); 

plot(z3,q3dall(l,:).','c','LineWidth',2); 
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plot(z3,Na3all(l,:).','b','LineWidth',2); 
plot(z3 ,C13 all( 1,:' ,'g', 'LineW idth' ,2); 

hold off; grid on; xlabel('z (nm)','FontSize',20); ylabel('Na, Cl, ArginineVFontSize',20); xlim([0 L+2*LR]); 
ylim([-0.1 ymaxl]); title('t=0 (a.u.)','FontSize',20); 
hh=legend('Na','Cl','c_l','c_2','c_3','c_4','Location',’SouthEasf); 
subplot(2,4,2); 

[tmp ,itonp3 ] =min(ab s (tv- (ton+3))); 

plot(zl,Nalall(itonp3,:).','b','LineWidth',2); hold on; 

plot(z 1 ,C11 all(itonp3,:) .','g','LineWidth',2); 

plot(z 1 ,q 1 aall(itonp3,:). ','r','LineWidth',2); 

plot(zl,qlball(itonp3,:).','k','LineWidth',2); 

plot(z 1 ,q 1 c all(itonp3,:).' ,'m','Line Width', 2); 

plot(zl,qldall(itonp3,:).','c','LineWidth',2); 

plot(z2,q2aall(itonp3,:).','r','LineWidth',2); 

plot(z2,q2ball(itonp3,:).','k','LineWidth',2); 

plot(z2,q2call(itonp3,:).','m','LineWidth',2); 

plot(z2,q2dall(itonp3,:).','c','LineWidth',2); 

plot(z3 ,q3aall(itonp3,:). ’,'r','LineWidth',2); 

plot(z3,q3ball(itonp3,:).','k','LineWidth',2); 

plot(z3,q3call(itonp3,:).','m','LineWidth',2); 

plot(z3,q3dall(itonp3,:).','c','LineWidth',2); 

plot(z3 ,Na3 all(itonp3,:).', 'b', 'LineW idth', 2); 

plot(z3,C13all(itonp3,:).','g','LineWidth',2); 

hold off; grid on; xlabel('z (nm)','FontSize',20); ylabel('Na, Cl, Arginine','FontSize',20); xlim([0 L+2*LR]); 

ylim([-0.1 ymaxl]); title(['t=' num2str(tv(itonp3))' (a.u.)'],''FontSize',20); 

subplot(2,4,3); 

[tmp,itoffm2]=min(abs(tv-(toff-2))); 

plot(zl,Nalall(itoffm2,:).','b','LineWidth',2); hold on; 

plot(z 1 ,C11 all(itoffm2,:).' ,'g', 'LineW idth' ,2); 

plot(zl,ql aall(itoffm2,:).' ,'r' ,'Line W idth', 2); 

plot(zl,q 1 ball(itoffm2,:).' ,'k', 'LineW idth', 2); 

plot(zl,qlc all(itoffm2,:).' ,'m', 'LineW idth', 2); 

plot(zl,qldall(itoffm2,:).','c','LineWidth',2); 

plot(z2,q2aall(itoffm2,:).','r','LineWidth',2); 

plot(z2,q2ball(itoffm2,:).','k','LineWidth',2); 

plot(z2,q2call(itoffm2,:).','m','LineWidth',2); 

plot(z2,q2dall(itoffm2,:).','c','LineWidth',2); 

plot(z3,q3aall(itoffm2,:).','r','LineWidth',2); 

plot(z3,q3ball(itoffm2,:).','k','LineWidth',2); 

plot(z3,q3call(itoffm2,:).','m','LineWidth',2); 

plot(z3,q3dall(itoffm2,:).','c','LineWidth',2); 

plot(z3,Na3all(itoffm2,:).','b','LineWidth',2); 

plot(z3,C13all(itoffm2,:).','g','LineWidth',2); 

hold off; grid on; xlabel('z (nm)','FontSize',20); ylabel('Na, Cl, Arginine','FontSize',20); xlim([0 L+2*LR]); 

ylim([-0.1 ymaxl]); title(['t=' num2str(tv(itoffm2))' (a.u.)'],'FontSize',20); 

subplot(2,4,4); 

plot(zl,Nalall(end,:).','b','LineWidth',2); hold on; 
plot(z 1 ,C11 all(end,:).','g','Line Width',2); 
plot(z 1 ,q 1 aall(end,:).' ,'r', 'Line Width', 2); 
plot(zl,q 1 ball(end,:).' ,'k', 'LineW idth', 2); 
plot(zl,qlc all(end,:).' ,'m', 'LineW idth', 2); 
plot(z 1 ,q ldall(end, :).','c','LineWidth',2); 
plot(z2 ,q2aall(end,:).' ,'r', 'LineW idth', 2); 
plot(z2,q2ball(end,:).','k','LineWidth',2); 
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plot(z2 ,q2c all(end,:).','m','Line Width', 2); 

plot(z2,q2dall(end,:).','c','LineWidth',2); 

plot(z3 ,q3 aall(end,:).' ,'r', 'Line W idth', 2); 

plot(z3,q3ball(end,:).','k','LineWidth',2); 

plot(z3,q3call(end,:).’,'m','LineWidth',2); 

plot(z3,q3dall(end,:).','c','LineWidth',2); 

plot(z3,Na3all(end,:).','b','LineWidth',2); 

plot(z3,C13all(end,:).','g','LineWidth',2); 

hold off; grid on; xlabel('z (nm)','FontSize',20); ylabel('Na, Cl, Arginine','FontSize',20); xlim([0 L+2*LR]); 

ylim([-0.1 ymaxl]); title(['t=' num2str(tv(end))' (a.u.)'],'FontSize',20); 

subplot(2,4,5); 

plot(zl,25*philall(l,:).','b','LineWidth',2); hold on; 

plot(z2,25 *phi2all( 1,:).', 'r' ,'LineWidth' ,2); 

plot(z3,25*phi3all(l,:).','b','LineWidth',2); hold off; 

grid on; xlabel('z (nm)','FontSize',20); ylabel('\phi (mV)','FontSize',20); 

xlim([0 L+2*LR]); ylim([(ymin2-0.5)*25 (ymax2+0.5)*25]); set(gca,'Ytick',(-4:l:0)*25); 

subplot(2,4,6); 

plot(zl,25*philall(itonp3,:).','b','LineWidth',2); hold on; 

plot(z2,25*phi2all(itonp3,:).','r','LineWidth',2); 

plot(z3,25*phi3all(itonp3,:).','b','LineWidth',2); hold off; 

grid on; xlabel('z (nm)','FontSize',20); ylabel('\phi (mV)','FontSize',20); 

xlim([0 L+2*LR]); ylim([(ymin2-0.5)*25 (ymax2+0.5)*25]); set(gca,'Ytick',(-4:l:0)*25); 

subplot(2,4,7); 

plot(zl,25*philall(itoffm2,:).','b','LineWidth',2); hold on; 

plot(z2,25*phi2all(itoffm2,:).','r','LineWidth',2); 

plot(z3,25*phi3all(itoffm2,:).','b','LineWidth',2); hold off; 

grid on; xlabel('z (nm)','FontSize',20); ylabel('\phi (mV)','FontSize',20); 

xlim([0 L+2*LR]); ylim([(ymin2-0.5)*25 (ymax2+0.5)*25]); set(gca,'Ytick',(-4:l:0)*25); 

subplot(2,4,8); 

plot(zl,25*philall(end,:).','b','LineWidth',2); hold on; 

plot(z2,25*phi2all(end,:).','r','LineWidth',2); 

plot(z3,25*phi3all(end,:).','b','LineWidth',2); hold off; 

grid on; xlabel('z (nm)','FontSize',20); ylabel('\phi (mV)','FontSize',20); 

xlim([0 L+2*LR]); ylim([(ymin2-0.5)*25 (ymax2+0.5)*25]); set(gca,'Ytick',(-4:l:0)*25); 

end 

iplot=input('Input 1 for showing animation, default is 1:'); 

if isempty(iplot) 

iplot=l; 

end 

if iplot==l 

int=input('Input animation time interval, default is 10: '); 

if isempty(int) 

int=10; 

end 

ipause=input('Input 1 to pause at each time frame, default is 0:'); 

if isempty(ipause) 

ipause=0; 

end 

for i=l:int:length(tv) 
figure(ll); 

subplot(2,2,l); plot(tv,phiLv,'Color','b','LineWidth',2); xlabel('f,'FontSize',18); ylabel('phiL','FontSize',18); grid on 
title(['t- num2str(tv(i))],'Fontsize',18); ylim([ymin2-l ymax2+l]); hold on; 
plot(tv(i),phiLv(i),'LineS tyle','None','Marker','.','MarkerSize',20,'Color','r'); 
set(gca,'Ytick',-6:2:6); hold off; xlim([0 tv(end)]); 
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subplot(2,2,2); plot(zLNalall(i,:).','b','LineWidth',2); hold on; 

plot(z 1 ,C11 all(i,:' ,'g', 'LineW idth' ,2); 

plot(z 1 ,q 1 aall(i, :).','r','LineWidth',2); 

plot(z 1 ,q 1 ball(i,:).','k','Line Width', 2); 

plot(z 1 ,q 1 c all(i,:).' ,'m', 'LineW idth' ,2); 

plot(zl,qldall(i,:).','c','LineWidth',2); 

plot(z2,q2aall(i,:).','r','LineWidth',2); 

plot(z2,q2ball(i,:).','k','LineWidth',2); 

plot(z2,q2call(i,:).','m','LineWidth',2); 

plot(z2,q2dall(i,:).','c','LineWidth',2); 

plot(z3,q3aall(i,:).','r','LineWidth',2); 

plot(z3,q3ball(i,:).','k','LineWidth',2); 

plot(z3 ,q3c all(i,:).'' ,'m', 'LineW idth', 2); 

plot(z3,q3dall(i,:).','c','LineWidth',2); 

plot(z3,Na3all(i,:).','b','LineWidth',2); 

plot(z3 ,C13 all(i,:).' ,'g', 'LineW idth' ,2); 

hold off; grid on; xlabel('z','FontSize',18); ylabel('Na, Cl, Arginine','FontSize',18); ylim([-0.1 l.l*ymaxl]); xlim([0 
L+2*LR]); 

% hh=legend('Na','Cl','Arg a'.'Arg b','Arg c'.'Aig d','Location','North'); 

% set(hh,'FontSize',8); 
subplot(2,2,3); 

plot(zl,philall(i,:).','k','LineWidth',2); hold on; 

plot(z2,phi2all(i,:).','k','LineWidth',2); 

plot(z3,phi3all(i,:).','k','LineWidth',2); hold off; 

grid on; xlabel('z','FontSize',18); ylabel('phi','FontSize',18); 

ylim([-6.5 6.5]); xlim([0 L+2*LR]); set(gca,'Ytick',-8:2:8); 

% set(gca,'Ytick',-60:20:60); 
if -isempty(ton) 
subplot(2,2,4) 

plot(tv,Imid,'Color','b','LineWidth',2); xlabel('t','FontSize',18); ylabel('Imiddle','FontSize',18); grid on; xlim([0 tv(end)]); 
% subplot(2,3,5) 

% plot(tv,Ql/V,'Color','b','LineWidth',2); xlabel('t','FontSize',18); ylabel('arginine fraction','FontSize',18); grid on; hold 
on; 

% plot(tv,Q3/V,'Color','g','LineWidth',2); plot(tv,Q2/V,'Color','r','LineWidth',2); hold off; 

% hh=legend('zone l','zone 3','zone 2'); set(hh,'FontSize',18); xlim([0 tv(end)]); 

% subplot(2,3,6) 

% plot(zl,Ileft(i)*ones(size(zl)),'b','LineWidth',2); xlabel('z','FontSize',18); ylabel('total current','FontSize',18); grid on; 
hold on; 

% plot(z2,Ileft(i)*ones(size(z2)),'g','LineWidth',2); plot(z3,Ileft(i)*ones(size(z3)),'r','LineWidth',2); 

% hold off; ylim([l.l*ymin3 l.l*ymax3]); xlim([0 L+2*LR]); 
end 

drawnow; 
if ipause==l 
pause; 
else 

pause(O.Ol); 

end 

end 

end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

function dudt=pnpld(t,u) 
global zl z2 z3 N1 N2 N3 Dzl Dz2 Dz3 
global DAzl DAz2 DAz3 Azl Az2 Az3 
global sphi K KS b S4offset LR L Qint 
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global NaL NaR C1L C1R phiL phiR phiLlow 
global DNa Dq Dq2 DC1 

global zNa zCl zq Lap Gammal Gamma2 Gamma3 
global iplot 

global ton toff Vwidtht 
global gij Wz 
Vwidth=Vwidtht; 
if isempty(ton) 
phiLactual=phiL; 
else 

phiLactual=((phiL-phiLlow)*tanh(Vwidth*(t-(ton+l)))-(phiL-phiLlow)*tanh(Vwidth*(t-(toff-l))))/2+phiLlow; 

% if t<ton I t>toff 
% phiLactual=phiLlow; 

% else 

% phiLactual=((phiL-phiLlow)*(l-exp(-Vwidth*(t-ton)))+(phiL-phiLlow)*(l-exp(-Vwidth*(toff-t))))/2+phiLlow; 

% end 
end 

Nal=u(l:Nl); Cll=u(Nl+l:2*Nl); qla=u(2*Nl+l:3*Nl); qlb=u(3*Nl+l:4*Nl); qlc=u(4*Nl+l:5*Nl); 
qld=u(5*Nl+l:6*Nl); stp=6*N 1; 

q2a=u(stp+l:stp+N2); q2b=u(stp+N2+l:stp+2*N2); q2c=u(stp+2*N2+l:stp+3*N2); q2d=u(stp+3*N2+l:stp+4*N2); 
stp=6*Nl+4*N2; 

Na3=u(stp+l:stp+N3); C13=u(stp+N3+l:stp+2*N3); q3a=u(stp+2*N3+l:stp+3*N3); q3b=u(stp+3*N3+l :stp+4*N3); 
q3c=u(stp+4*N3+l:stp+5*N3); q3d=u(stp+5*N3+l:stp+6*N3); S4disp=u(end); 

Nazl=Dzl*Nal; Naz3=Dz3*Na3; Clzl=Dzl*Cll; Clz3=Dz3*C13; 
qzla=Dzl*qla; qzlb=Dzl*qlb; qzlc=Dzl*qlc; qzld=Dzl*qld; 
qz2a=Dz2*q2a; qz2b=Dz2*q2b; qz2c=Dz2*q2c; qz2d=Dz2*q2d; 
qz3a=Dz3*q3a; qz3b=Dz3*q3b; qz3c=Dz3*q3c; qz3d=Dz3*q3d; 

% electric field: 

rhs 1 =-(zNa*Na 1 +zCl*Cl 1 +zq*(q 1 a+q 1 b+q 1 c+q 1 d))/Gamma 1; rhsl(l)=phiLactual; rhsl(end)=0; 
rhs2=-(zq*(q2a+q2b+q2c+q2d))/Gamma2; rhs2([l end])=0; 

rhs3=-(zNa*Na3+zCl*C13+zq*(q3a+q3b+q3c+q3d))/Gamma3; rhs3(l)=0; rhs3(end)=phiR; 
rhs=[rhsl;rhs2;rhs3]; warning on; phi=Lap\rhs; warning off; 
phil=phi(l:Nl); phi2=phi(Nl+l:Nl+N2); phi3=phi(Nl+N2+l:end); 
phizl=Dzl*phil; phiz2=Dz2*phi2; phiz3=Dz3*phi3; 

% trap: 
warning off; 

tmp=Azl(l)*qla.*zl; tmp=Dzl\tmp; tmpl=tmp(end)-tmp(l); 
tmp=Az2(l)*q2a.*z2; tmp=Dz2\tmp; tmp2=tmp(end)-tmp(l); 
tmp=Az3(l)*q3a.*z3; tmp=Dz3\tmp; tmp3=tmp(end)-tmp(l); 
qapos=(tmp 1 +tmp2+tmp3)/Qint; 

tmp=Azl(l)*qlb.*zl; tmp=Dzl\tmp; tmpl=tmp(end)-tmp(l); 
tmp=Az2(l)*q2b.*z2; tmp=Dz2\tmp; tmp2=tmp(end)-tmp(l); 
tmp=Az3(l)*q3b.*z3; tmp=Dz3\tmp; tmp3=tmp(end)-tmp(l); 
qbpos=(tmpl+tmp2+tmp3)/Qint; 

tmp=Azl(l)*qlc.*zl; tmp=Dzl\tmp; tmpl=tmp(end)-tmp(l); 
tmp=Az2(l)*q2c.*z2; tmp=Dz2\tmp; tmp2=tmp(end)-tmp(l); 
tmp=Az3(l)*q3c.*z3; tmp=Dz3\tmp; tmp3=tmp(end)-tmp(l); 
qcpos=(tmp 1 +tmp2+tmp3)/Qint; 

tmp=Azl(l)*qld.*zl; tmp=Dzl\tmp; tmpl=tmp(end)-tmp(l); 
tmp=Az2(l)*q2d.*z2; tmp=Dz2\tmp; tmp2=tmp(end)-tmp(l); 
tmp=Az3(l)*q3d.*z3; tmp=Dz3\tmp; tmp3=tmp(end)-tmp(l); 
qdpos=(tmpl+tmp2+tmp3)/Qint; 
warning on; 

qarel=qapos-(S4disp-3*sphi); qbrel=qbpos-(S4disp-sphi); qcrel=qcpos-(S4disp+sphi); qdrel=qdpos-(S4disp+3*sphi); 
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dS4dispdt=(-KS*(S4disp-(LR+L/2+S4offset))+K*(qarel+qbrel+qcrel+qdrel))/b; 
Vzla=K*(zl-(S4disp-3*sphi)); Vzlb=K*(zl-(S4disp-sphi)); Vzlc=K*(zl-(S4disp+sphi)); Vzld=K*(zl 
(S4disp+3*sphi)); 

Vz2a=K*(z2-(S4disp-3*sphi)); Vz2b=K*(z2-(S4disp-sphi)); Vz2c=K*(z2-(S4disp+sphi)); Vz2d=K*(z2 
(S4disp+3*sphi)); 

Vz3a=K*(z3-(S4disp-3*sphi)); Vz3b=K*(z3-(S4disp-sphi)); Vz3c=K*(z3-(S4disp+sphi)); Vz3d=K*(z3 
(S4disp+3*sphi)); 

% flux: 

JNazl=-Azl.*DNa.*(Nazl+zNa*Nal.*phizl); 

JClz 1=-Az 1. *DC1. *(Clz 1 +zCl*Cl 1. *phiz 1); 

Jqzla=-Azl.*Dq.*(qzla+zq*qla.*phizl+qla.*Vzla+gij*qla.*(qzlb+qzlc+qzld)); 

Jqzlb=-Azl.*Dq.*(qzlb+zq*qlb.*phizl+qlb.*Vzlb+gij*qlb.*(qzla+qzlc+qzld)); 

Jqzlc=-Azl.*Dq.*(qzlc+zq*qlc.*phizl+qlc.*Vzlc+gij*qlc.*(qzla+qzlb+qzld)); 

Jqz 1 d=-Az 1. *Dq. *(qz 1 d+zq*q 1 d. *phiz 1 +q 1 d. * Vz 1 d+gij *q 1 d. *(qz 1 a+qz 1 b+qz 1 c)); 

Jqz2a=-Az2.*Dq2.*(qz2a+zq*q2a.*phiz2+q2a.*(Wz+Vz2a)+gij*q2a.*(qz2b+qz2c+qz2d)); 

Jqz2b=-Az2.*Dq2.*(qz2b+zq*q2b.*phiz2+q2b.*(Wz+Vz2b)+gij*q2b.*(qz2a+qz2c+qz2d)); 

Jqz2c=-Az2.*Dq2.*(qz2c+zq*q2c.*phiz2+q2c.*(Wz+Vz2c)+gij*q2c.*(qz2a+qz2b+qz2d)); 

Jqz2d=-Az2.*Dq2.*(qz2d+zq*q2d.*phiz2+q2d.*(Wz+Vz2d)+gij*q2d.*(qz2a+qz2b+qz2c)); 

JNaz3=-Az3.*DNa.*(Naz3+zNa*Na3.*phiz3); 

JClz3=-Az3.*DC1.*(Clz3+zCl*C13 *phiz3); 

Jqz3a=-Az3.*Dq.*(qz3a+zq*q3a.*phiz3+q3a.*Vz3a+gij*q3a.*(qz3b+qz3c+qz3d)); 

Jqz3b=-Az3.*Dq.*(qz3b+zq*q3b.*phiz3+q3b.*Vz3b+gij*q3b.*(qz3a+qz3c+qz3d)); 

Jqz3c=-Az3.*Dq.*(qz3c+zq*q3c.*phiz3+q3c.*Vz3c+gij*q3c.*(qz3a+qz3b+qz3d)); 

Jqz3d=-Az3.*Dq.*(qz3d+zq*q3d.*phiz3+q3d.*Vz3d+gij*q3d.*(qz3a+qz3b+qz3c)); 

% conservation law: 

dNadtl=-DAzl*JNazl; dNadt3=-DAz3*JNaz3; dCldtl=-DAzl*JClzl; dCldt3=-DAz3*JClz3; 
dqdtla=-DAzl*Jqzla; dqdtlb=-DAzl*Jqzlb; dqdtlc=-DAzl*Jqzlc; dqdtld=-DAzl*Jqzld; 
dqdt2a=-DAz2*Jqz2a; dqdt2b=-DAz2*Jqz2b; dqdt2c=-DAz2*Jqz2c; dqdt2d=-DAz2*Jqz2d; 
dqdt3a=-DAz3*Jqz3a; dqdt3b=-DAz3*Jqz3b; dqdt3c=-DAz3*Jqz3c; dqdt3d=-DAz3*Jqz3d; 

% BC: 

dNadtl(l)=Nal(l)-NaL; 
dN adt 1 (end)=JN az 1 (end); 
dNadt3(l)=JNaz3(l); 
dN adt3 (end)=N a3 (end) -NaR; 
dCldt 1(1 )=C11(1 )-ClL; 
dCldt 1 (end)=JClz 1 (end); 
dCldt3 (l)=JClz3(l); 
dCldt3(end)=C13(end)-ClR; 

dqdt 1 a( 1 )=Jqz 1 a( 1); dqdt 1 a(end)=q2a( 1 )-q 1 a(end); 

dqdt2a( 1 )=Jqz2a( 1)-Jqz 1 a(end); dqdt2a(end)=q2a(end)-q3a( 1); 
dqdt3a( 1 J=Jqz3a( 1 )-Jqz2a(end); dqdt3a(end)=Jqz3a(end); 
dqdt 1 b( 1 )=Jqz lb( 1); dqdt 1 b(end)=q2b( 1 )-q 1 b(end); 

dqdt2b( 1 )=Jqz2b( 1)-Jqzlb(end); dqdt2b(end)=q2b(end)-q3b( 1); 
dqdt3b(l)=Jqz3b(l)-Jqz2b(end); dqdt3b(end)=Jqz3b(end); 
dqdt lc(l)=Jqzlc(l); dqdt 1 c(end)=q2c( 1 )-q 1 c(end); 

dqdt2c( 1 )=Jqz2c( 1)-Jqzlc(end); dqdt2c(end)=q2c(end)-q3c( 1); 
dqdt3c( 1 )=Jqz3c( 1)-Jqz2c(end); dqdt3c(end)=Jqz3c(end); 
dqdt 1 d( 1 )=Jqz 1 d( 1); dqdt 1 d(end)=q2d( 1 )-q 1 d(end); 

dqdt2d( 1 )=Jqz2d( 1)-Jqzld(end); dqdt2d(end)=q2d(end)-q3d( 1); 
dqdt3d(l)=Jqz3d(l)-Jqz2d(end); dqdt3d(end)=Jqz3d(end); 

dudt=[dNadt 1 ;dCldt 1 ;dqdt 1 a;dqdt 1 b ;dqdt 1 c ;dqdt 1 d;dqdt2a;dqdt2b ;dqdt2c ;dqdt2d; ... 

dNadt3;dCldt3;dqdt3a;dqdt3b;dqdt3c;dqdt3d;dS4dispdt]; 

[t max(abs(dudt))] 
if iplot==l & t>0 
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figure(50); subplot(l,2,l); 

plot(zl,Nal,'b','LineWidth',2); hold on; plot(z3,Na3,'b','LineWidth',2); plot(zl,Cll,'g','LineWidth',2); 
plot(z3,C13,'g','LineWidth',2); 

plot(zl,qla,'rVLineWidth',2); plot(zl,qlb,Y,TdneWidth',2); plot(zl,qlc,'mVLineWidth',2); 
plot(zl,qld,'c','LineWidth',2); 

plot(z2,q2a,'r','LineWidth',2); plot(z2,q2b,'kVLineWidth',2); plot(z2,q2c,'m','LineWidth',2); 
plot(z2,q2d,'c','LineWidth',2); 

plot(z3,q3a,'r','LineWidth',2); plot(z3,q3b,'k','LineWidth',2); plot(z3,q3c,'m','LineWidth',2); 
plot(z3,q3d,'c','LineWidth',2); 

hold off; grid on; xlabel('z'); ylabel('Na, Cl, q'); title(['t=' num2str(t)]); 
subplot( 1,2,2); 

plot(zl,phil,'k','LineWidth',2); hold on; plot(z2,phi2,'k','LineWidth',2); plot(z3,phi3,'k','LineWidth',2); 

grid on; xlabel('z'); ylabel('\phi'); hold off; drawnow; 

end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% CHEB compute D = differentiation matrix, x = Chebyshev grid 

function [D,x] = cheb(N) 
if N==0, D=0; x=l; return, end 
x = cos(pi*(0:N)/N)'; 
c = [2; ones(N-l,l); 2].*(-l). A (0:N)'; 

X = repmat(x,l,N+l); 
dX = X-X'; 

D = (c*(l./c)')./(dX+(eye(N+l))); % off-diagonal entries 

D = D - diag(sum(D')); % diagonal entries 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%% 

function [x, DM] = chebdif(N, M) 

% The function [x, DM] = chebdif(N,M) computes the differentiation 
% matrices Dl, D2,..., DM on Chebyshev nodes. 

% 

% Input: 

% N: Size of differentiation matrix. 

% M: Number of derivatives required (integer). 

% Note: 0<M<=N-1. 

% 

% Output: 

% DM: DM(l:N,l:N,ell) contains ell-th derivative matrix, ell=l..M. 

% 

% The code implements two strategies for enhanced 
% accuracy suggested by W. Don and S. Solomonoff in 
% SIAM J. Sci. Comp. Vol. 6, pp. 1253-1268 (1994). 

% The two strategies are (a) the use of trigonometric 
% identities to avoid the computation of differences 
% x(k)-x(j) and (b) the use of the "flipping trick" 

% which is necessary since sin t can be computed to high 
% relative precision when t is small whereas sin (pi-t) cannot. 

% Note added May 2003: It may, in fact, be slightly better not to 
% implement the strategies (a) and (b). Please consult the following 
% paper for details: "Spectral Differencing with a Twist", by 
% R. Baltensperger and M.R. Trummer, to appear in SIAM J. Sci. Comp. 

% J.A.C. Weideman, S.C. Reddy 1998. Help notes modified by 
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% JACW, May 2003. 


I = eye(N); % Identity matrix. 

L = logical(I); % Logical identity matrix. 

nl = floor(N/2); n2 = ceil(N/2); % Indices used for flipping trick. 

k = [0:N-1]'; % Compute theta vector, 

th = k*pi/(N-l); 


x = sin(pi*[N-l :-2:1 -NJ'/(2*(N-1))); % Compute Chebyshev points. 

T = repmat(th/2,l,N); 

DX = 2*sin(T'+T).*sin(T'-T); % Trigonometric identity. 

DX = [DX(l:nl,:); -flipud(fliplr(DX(l:n2,:)))]; % Flipping trick. 

DX(L) = ones(N,l); % Put l's on the main diagonal of DX. 

C = toeplitz((-l). A k); % C is the matrix with 

C(l,:) = C(l,:)*2; C(N,:) = C(N,:)*2; % entries c(k)/c(j) 

C(:,l) = C(:,l)/2; C(:,N) = C(:,N)/2; 

Z = l./DX; % Z contains entries l/(x(k)-x(j)) 

Z(L) = zeros(N,l); % with zeros on the diagonal. 

D = eye(N); % D contains diff. matrices. 

for ell = 1:M 

D = ell*Z.*(C.*repmat(diag(D),l,N) - D); % Off-diagonals 
D(L) = -sum(D'); % Correct main diagonal of D 

DM(:,:,ell) = D; % Store current D in DM 

end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

function F = myfun2(x,xdata) 

F = x(l)*exp(-xdata/x(2))+x(3)*exp(-xdata/x(4)); 
function F = myfun3(x,xdata) 

F = x(l)*exp(-xdata/x(2))+x(3)*exp(-xdata/x(4))+x(5)*exp(-xdata/x(6)); 
function F = myfun4(x,xdata) 

F = x(l)*exp(-xdata/x(2))+x(3)*exp(-xdata/x(4))-i-x(5)*exp(-xdata/x(6))-i-x(7)*exp(-xdata/x(8)); 
function F = myfun5(x,xdata) 

F = x(l)*xdata. A x(5).*exp(-xdata/x(2))+x(3)*xdata. A x(6).*exp(-xdata/x(4)); 
function F = myfun6(x,xdata) 

F = x(l)*xdata. A x(7).*exp(-xdata/x(2))+x(3)*xdata. A x(8).*exp(-xdata/x(4))+x(5)*xdata. A x(9).*exp(-xdata/x(6)); 
function F = myfun7(x,xdata) 

F = x(l)*xdata. A x(9).*exp(-xdata/x(2))+x(3)*xdata. A x(10).*exp(-xdata/x(4))-i-x(5)*xdata. A x(ll).*exp(- 
xdata/x(6))+x(7)*xdata. A x(12).*exp(-xdata/x(8)); 
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